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Inversion of Symmetric Coefficient Matrix of 
Positive-Definite Quadratic Form 


1. Introduction. This paper presents a method of inverting a positive-definite 
symmetric matrix which has been coded by the author for the IBM 704 digital 
computer at Los Alamos. The code is now running so that timing and accuracy 
descriptions are included. 

For brevity we shall call a real » by m symmetric matrix a positive-definite if 
the quadratic form x7ax is positive-definite, x real. We shall invert a by construct- 
ing an upper triangular matrix H subject to the conditions |H|=1 and 
H'aH = d—, where d~ is diagonal. It then follows that a! = HdH™, |a| =|d“|, 
and x"’ax = ¢7d“£, where x = Hé. This gives the 


Criterion: The necessary and sufficient conditions that x7ax be positive- 
definite are that alJ » terms of the diagonal matrix d-', which the code forms for 
printing, be positive. 

It should be remarked that, while the positive-definite condition ensures no 
attempted machine division by zero, there exists a large class of indefinite matrices 
which can be inverted by our method. 

2. The choice of H. The method of constructing H so as to diagonalize a 
according to H™aH = d- will be described for the case m = 4, the generalization 
to any » being immediate. We employ Gauss elimination to define 


(2.1) hi Ss - @1;/411, ay = (Gudp = Gy4;;)/a1, < 2 = 2, 3, 4, 


the division being justified by our positive-definite hypothesis on a, and form 


1 hin fig hus ay 0 0 0 


A, = : ; : : | a® = H,'aH, = 


00 0 1 


2) 2) 
0 @22 ays aX 
0 oP of of 
0 


(2) 


a a? af 


for which |H;| = 1. Since x"ax = y™a®y, where x = Hyy, the right-hand form 
is likewise positive-definite. 
Continuing with the diagonalization by defining 


LS (2) 7 (2) 3) _ /,(2),(2) (2) (2), / (2) 
ha; = — a3j/ay, af’ = fay’ — agay’)/ay’, 


for j, k = 3, 4, and so on, we finally obtain 
his hus 


ooor 
ooro 
coor 
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Iye-1 era-hees + Ieig-1 era (eoshtsa + eas) + Aas -ha + Aria-1| 
e 1 hes-1 hos hag + haat 
(2.2) H= 0 i heat 
0 0 1 


d— = diag. matrix of elements (a1:, a3, a§}, a). 

An inspection of the inner product composition of the columns of H shows 
how this matrix generalizes for arbitrary . H is generated from the triangular 
array of the h’s by first forming the elements of the last column, proceeding from 
the bottom up, then passing to the adjacent column, etc. 

It may be helpful to point out the relation which the method of this section 
bears to a familiar procedure for reducing a positive-definite quadratic form x7ax 
to a sum of squares. Namely, one employs the algorithm of Gauss-Banachiewicz 
(see Bodewig [2], p. 105, or Turing [3], p. 289), to factor the coefficient matrix 
according to the form 


(2.3) a= (I+ U)*D(I + U), 


where (I + U) is a uniquely determined upper triangular matrix with diagonal 
elements unity and D is a diagonal matrix of positive terms. Comparison with 
our factorization H'aH = d- identifies H with (J + U)— and d— with D. Then 
the transformation H-'x = £ gives the desired reduction to a sum of squares, 
x?ax = §™D¢. Thus, if we had chosen to employ the algorithm (2.3) to form 
H- = (I + U), we would then have been confronted, as was Banachiewicz [4], 
p. 398, with the additional task of inverting the transformation H-'x = é to 
express it in the solved form x = H¢é necessary for our matrix inversion scheme. 

It can be shown [2], p. 91, that the elements d;—' of the diagonal matrix d- 
are given by the quotients of successive principal minors of a, 


dy" =|au|, det =|ai022|/|au|, ds = |@1022033|/|@nde], ---. 


3. The inverse matrix a. The matrix H of (2.2) is of the upper triangular 
form 


M11 «712 8 714 
0 

(3.1) H= 0 “Y or > M1 = 22 = 233 = na = 1, 
© "44 


so that the distinct elements of the symmetric inverse matrix a~! = HdH” are 
given by 


aden Medre Mdm nde 

3.2 —l = Neder Naedese Noedenes 
. nscdnse nsedanas 
nad anes 


r=1,2,3,4 s=2,3,4 ¢=3,4. 
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In writing (3.2) the convention of summing the repeated indices r, s, t, over their 
respective ranges has been observed. 
4. Storage arrangement for coding. The input array is illustrated for » = 3: 


Go(@11) Gi@i2) G2(a:s) 
Gs(@22) Ga(ass) (3n(" + 1) words) 
(4.1) Gs(@ss) 
Ge(b1) Gr(b2) Gs(bs) (n words) 
IGo( ) Gl ) Gul ) Gil )|| (+ 1 words reserved for output) 


where the G’s represent $n(m + 1) + 2n+ 1 consecutive storage locations 
The b's are the right-hand sides of any linear system ax = 5, while the contents 
of the last » + 1 words at the time of input are immaterial. 

A code of 263 orders has been written which requires a fixed block of 13 and a 
variable block of 2” erasable words in the consecutive locations Dg, - --, Dan—1. 
Since the computations of the type (2.1) involve only two rows at a time, these 
two rows can be placed in the D-block for ease in addressing and the results of the 
computation stored appropriately in the G-block, destroying data no longer 
required. 

Proceeding in this way, the triangular array (3.1) for H is generated in the 
first 4n(n + 1) locations of the G-block, the former non-unity diagonal terms 
which define d-' having been stored for printing and further use in the first ” 
locations of the last line of (4.1). 

Inspection of the matrix (3.2) defining a~ reveals that again only one pair of 
rows at a time is involved in forming a". Once more the D-block serves to store 
this pair for easy addressing and clearance in storing a in the first $n(m + 1) 
locations of the G-block, data no longer needed being destroyed. 

The result is the output array, illustrated for » = 3, available for printing or 
further computing : 


Golan) Gilerw)  Ge(ars) 
Gs(a22)  Ge(aes) (inverse matrix a~) 
(4.2) Gs(ass) 
Ge(x1)  Gr(x2) Gs (xs) (solution of ax = 5) 
Ge(di") Gr(ds") Guilds) Gis(\a|)|| (pos.-def. criterion and det.). 


5. Computing with a~. It is anticipated that the elements a; of a~' may be 
subsequently employed to form inner products of the type = ajryr, where y is 


an arbitrary vector. To facilitate the addressing of this cum, ion may load the 
jth row of a into Do, ---, D,-1 and the y’s into D,, ---, Deas. The inversion 
code contains a block of 20 consecutive orders which iether the jth row of a 
from the triangular output’array of (4.2) and store in Do, ---, D,-1. After the 
inversion these 20 orders may be retained in memory for this purpose. 

6. Accuracy, timing and limitations on nm. The Hilbert matrix H,, 

= (“¢+ 7 — 1)", i, 7 = 1, 2, ---, m, being somewhat a champion among ill- 
conditioned matrices, imposes a severe test upon the accuracy of any method of 
inversion. The true inverse of H, has been tabulated by I. R. Savage and E. Lukacs 





58 INVERSION OF SYMMETRIC COEFFICIENT MATRIX 


({1], p. 105), through » = 10. A floating-point nine-digit approximation to H,, 
which took account of all digits in the IBM 704, was loaded and the computed 
inverse of this approximate H, was compared with the true inverse of H,. The 
following table of comparisons was obtained : 


| H,| Agreement with true H,- 


1.6534(10-") 5 digits 

3.75(10-”) 3 digits 

5.3(10-") 2 digits 

4(10-*5) max. discrepancy 20% 

2(10-* max. discrepancy 50% 
10-* inversion failed. 


Conraur B 


J. Todd reported in 1954 ([1], p. 113), that SEAC failed to invert H, when using 
a single precision fixed point routine. 

A less stringent test matrix isT,, yi; = yin = —t(n +1 — 7)/(n + 1),7 <j, 
whose inverse ([1], p. 112), iscy = — 2,¢3 = 1, |¢ —j|=1,¢,; = 0, |¢ — j|>1, 
t,j = 1,2, ---,”. A floating-point, eight-digit approximation to I, was loaded 
and the computed inverse of this approximate T,, was compared with the true 
T,,'. The following two runs were made: 


Agreement with 
n IT.| true T,— Computing time 
49 — 0.020000 6 digits 1min 33 sec 
115 — 0.00862 5 digits 19 min 30 sec. 


The computing time for inversion seems to be proportional to m*. Comparison 
with SEAC’s time of about three hours to invert I'ys, as reported by J. Todd in 
1954 [1], p. 113, gives an idea of the progress being made in the design of 
high-speed computers. 
For a machine with 2° words of core storage, the inequality which limits the 
_size of n is $n(n + 9) + 1 < 2° — 295, where 295 is the storage required by the 
Los Alamos print program. This gives the table: 


2¢ Max. n Approx. running time 


4,096 82 7 min 
8,192 121 23 min 
32,768 250 200 min 


Tuomas C. DOYLE 


University of California 
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NUMERICAL EVALUATION OF MULTIPLE INTEGRALS I 


Numerical Evaluation of Multiple Integrals I 


Introduction. Several specific methods for numerical evaluation of integrals 
over higher dimensional regions have been proposed. James Clerk-Maxwell [1] 
proposed the formulas for the rectangle and the rectangular parallelopipedon in 
1877. Appell, Burnside, Ionescue, and Mineur have developed special formulas 
for planar regions. Tyler [2] recently gave formulas for rectangles, parabolic 
regions, cubes, and for the hypercubes. Others have developed formulas pri- 
marily for rectangular regions based on the formulas for the line. 

Although the original draft of this paper antedates their results, Hammer, 
Marlowe and Stroud [3] have given an inductive method for numerical evaluation 
of integrals precise for kth degree polynomials over the n-simplex. They also show 
how to obtain formulas for cones based on regions for which integration formulas 
are given. Certain affinely symmetrical formulas for triangle and tetrahedron are 
given in their paper. In a sequel Hammer and Stroud [4] gave affinely sym- 
metrical formulas precise for the quadratic and cubic polynomials over the 
n-simplex. W. H. Peirce [5] in his doctoral dissertation has given numerical 
integration formulas over the circular annulus precise for arbitrarily high degree 
polynomials. 

Despite the variety of results and their particular interest, the theory of 
numerical integration in higher dimensions is in a very crude state of development 
in comparison with numerical integration of functions of one variable. For this 
reason we give here a few theorems which are actually quite simple but which 
evidently have escaped the awareness of most research workers. These theorems 
form a partial basis for development of classical type integration formulas. 

The sole numerical integration device with claims to somewhat universal 
applicability is that called the Monte Carlo method which has been extended 
and popularized by S. M. Ulam and John von Neumann. Granted the presence 
of a suitable sequence of ‘‘random”’ digits or the means of generating such, the 
Monte Carlo methods can be applied to a large variety of problems by using 
high-speed digital computers. It has been stated as an advantage of Monte Carlo 
methods that the number of evaluation points for, say, 5 percent error does not 
mount as rapidly for higher dimensions as it does for classical formulas. We shall 
show (as Tyler, and Hammer and Stroud have shown in other cases) that the 
evaluation points for classical type formulas can conceivably be kept down very 
well. However, the problem of determining points and weights in these circum- 
stances is a problem in the theory of equations which is difficult inherently and 
also because of the great variety of functions and regions of conceivable interest. 
In this direction the Monte Carlo method is less sensitive to the deviations among 
regions. 

A Basic Theorem for Numerical Integration. In this section we state a theorem 
which greatly extends the usefulness of particular integration formulas when they 
are available. We limit the discussion to integration formulas of the form 


(1) f w(x)f(xde = F arf(x), 
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where R is an n-dimensional region in Euclidean n-space, E,, x is a point vector in 
E,, 41, ---,@m are numbers (real or complex) and evaluation points x1, ---, %» 
are in the domain of f. In particular, R will be taken as the closure of a bounded 
open set in Z, and the weight function w(x) and f are continuous on R. That these 
conditions may be relaxed will not concern us. The dx indicates a volume ele- 
ment in £,. 

For each relevant f and a given formula (1) there is the associated error, a 
number E(f), which we define as 


(2) E(f) = Daf) - f w(x) f(x)de. 


Let S be another region in EZ, such that there exists a transformation J with 
continuous nonvanishing Jacobian which maps R onto .S. We write y = Jx for 
x € R. Then if g(y) is defined and continuous on S we may write 


~ Pes oT 


(3) ff woreordy = f w(y)g(9) | Jae, 


‘SITES EBes 


where in the integrand on the right we have y = Jx and || is the absolute value 
of the Jacobian which may also be indicated as dy/dx. Now if we have an inte- 
gration formula (1) for R and the choice w(x) = w:(y), then with f(x) replaced 
in (2) by |J|g(y) we have 


one coo 


(4) E(|Jlg) = XaslJileos) — f we0)| lax, 


where y; = Jx; and |J;| is the value of |J| at x;. ~ } 


THEOREM 1. Writing 


(5) E(S, 2) = Xbay,) — f w,(y)e(y)dy, 


where b; = a;|J;|, then E(S, g) = E(|J|g) and hence from every formula of type 
(1) for R there is a corresponding formula of the same type for S with error func- 
tienal E(S, g). 








Coro.iary 1.1. If |J| is constant on R then E(S, g) =|J\|E(g) and hence if 
E(g) vanishes for a function g so does E(S, g). In particular, if the transformation J 
is a nonsingular affine transformation, |J| is constant. 

CorOLLary 1.2. Identifying w(x) with w;(y)|J|, y = Jx and f(x) with g(y), 
assuming a formula (1), gives from equation (2) 







(6) E(S, 2) = E(g) = Deg(y,) - f w:(y)e(y)dy. 





Theorem 1 and the two corollaries make it possible to generate integration 
formulas for a variety of regions. The simplest and the most important case arises 
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when J is an affine transformation (Corollary 1.1). Thus certain numerical inte- 
gration formulas over a particular sphere precise for polynomials of at most 
degree k give immediately integration formulas for any ellipsoid precise for 
polynomials of at most degree k. 

Cartesian Product Regions. Let R be the Cartesian product of two regions 
R, and R; in lower dimensional Euclidean spaces. Let us designate y « Ri, z¢« Rz 
so that every x « R can be written (y, 2). Suppose there are numerical integration 
formulas given of type (1) for weight functions w;(y) and we(z), defined on R, 
and R, respectively, so that the error functionals may be written 


(1) E(R: fr) = Dafiy) — f ws (y)fily)dy 


(2) E(Rs, fs) = X bifale)) — f ws(2)fo(2)ds. 


Then we ask for a numerical integration formula over R = R, X R2z with weight 
function w:(y)we(z). We have, by using (1) and (2), 


THEOREM 2. Let f(x) = f(y, 2) be a function defined and continuous on R = R, 
X Re. Then if we define 


3) ERA) =CCedsns) — f f wi (y)we(s) f(y, #)dyds, 


we have 


E(R, f) = 5 bELRy f(y, 5))] + f w:(y)ELRs, f(y, #) dy 


= 5 aELRs, flys, 2)] + f w,(2)ELR,, f(y, 2) Ws. 


CorOLiary 2.1. If f(y, 2) is a function such that ELR,, f(y, 2)] = 0 for each 
ze R; and E[R:, f(y, z)] = 0 for each y « R,; then E(R, f) = 0 and the formula (3) 
of Theorem 2 is exact for f. 

If F; is a class of functions defined on Ri, and F, is a class of functions defined 
on Ro, then the Cartesian product class, F = F,; X F2, of functions defined on 
Ri X R: is the class of all functions f(y, z) such that f(y, 2) « F; for each ze Re 
and f(y, 2) « F2 for each y e Ri. 

Coro.uary 2.2. If E(Ri, f:) = 0 for fie Fi: and E(Re, fe) = 0 for foe F2 then 
necessarily E(R, f) = 0 for every f « Fi X F2, where the error functionals are defined 
by (1), (2), (3). Im other words, the numerical integration formula in (3) is neces- 
sarily exact for the Cartesian product of the classes of functions for which (1) and 
(2) are exact. 

CoroLiary 2.3. If classes of functions F and G respectively defined on R, and 
R; are representable as all linear combinations of basis sets of functions fi, ---, fp 
and gi, ---, g, respectively then F X G is the set of all functions with (fig;) as basis. 

Proof: lfh = > > cifig; where the c;; are constants, then he F X G. To show 


‘ ? 
that he F X G implies fA is a linear combination of the fig;, we observe that 
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h = > asf; = > b,g; where the a; are unique functions on R:, and b; are unique 


s 2 
functions on R,. Let 1, y2, ---, ¥p be the points in R,; such that the determinant 
| f:(yx) | is not zero. Then we have with bj the value of 5; at y: 


L afi(ye) ms L bagi k=1,---,p 


which are identities on R». Solving for ai, ---,@, we find a; = >> cig; where ¢;; 
7 
are constants and are unique. Hence h = > > cijfig;. 
ii 


' The usefulness of the Cartesian product principle is extensive. However, even 
here there must be limitations. Thus over the hypercube in E, the generation of 
a formula based on Cartesian powers of the line segment formulas with m points 
gives m” points. If ” is large and m > 2, this will involve a large number of points 
of evaluation. 


Examples. Gauss’ four-point formula on the line is exact for at most 7th degree 
polynomials in one variable. Hence the 16-point formula 


EE aah, x) = f f #00, »)aedy 
R 


where a; are the Gauss weights, R is the square, and x; are the zeros of the Le- 
gendre polynomial, is exact for f(x,y) = > Lexy? O0<i<7, O<j7 <7 
since this class of polynomials is the Cartesian square of the class of seventh 
degree polynomials in one variable. Observe that this class contains all seventh 
degree polynomials and additional polynomials since 7 and 7 range independently 
from 0 to 7. Similarly formulas may be developed for the hypercubes. 

Laguerre’s 3-point formula for weight function e~*, interval 0, © holds for 
all fifth degree polynomials. Hence 


3 3 oO C) 
Easley x) = f f eet f(x, y)dedy 
1 1 0 0 


holds exactly for every f of form >> > c;;x*y’. Thus we obtain a formula for the 
iS<5 755 

quadrant. A linear transformation will give a formula for an infinite planar sector. 

The x; are the zeros of the Laguerre polynomials. Again these may be used to 

obtain formulas for octants, etc., in higher dimensions. 

Combinations of Gauss formulas and Laguerre formulas give formulas for 
the planar half-strip, and combinations of Hermite’s and Laguerre’s will give 
formulas for the half-plane. A formula for a planar region generates a formula 
for the cylinder based on this region. In scattering and potential problems Car- 
tesian product regions which are powers or products of 3-dimensional regions are 
common. Thus formulas for each region are useful for generating others. 

P. Davis of the National Bureau of Standards wrote of an example he and 
P. Rabinowitz estimated and which has since been published [6 ]. This example is 
to integrate exp (x1%2%3x4) over the unit cube in 4-space. We calculated an esti- 
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mate using the fourth ‘“‘power’’ of the Gauss two-point formula, and obtained 
with these 16 points half the error reported by Davis and Rabinowitz using 
2'* = 32,768 points and a Monte Carlo method. Since the cost is roughly propor- 
tional to the number of points, it is seen that the classical type formula is here 
much more efficient. Of course this function is particularly fortunate for Cartesian 
power methods since it is expandable in powers of x:xoxsx4. 

Particular Formulas for Symmetrical Regions. In a sequel to this paper we 
will publish a table of formulas for the regions we have considered, in detail. 
However, we prefer to define symmetrical regions and correspondingly sym- 
metrical formulas in this paper and to illustrate applications in one type of case 
since consideration of these problems has led to a variety of solved and ‘unsolved 
related ones. 

Rather than introduce specialized terms we shall say a region R in n-space is 
symmetrical if : whenever it contains a point x, it also contains all points obtained 
from x by interchanging coordinates and changing signs of coordinates. Thus in 
general with a point x there are required 2* m! points (including x). A symmetrical 
region R is invariant under linear transformations with matrix containing exactly 
one non-zero element in each row and column, that non-zero element being either 
+1 or —1. In evaluating integrals of monomials over R, symmetry as defined, 
saves effort as indicated in the following theorem. 


THEOREM 3. The integral over a symmetrical region R of any monomial containing 
an odd power is zero, the integral of a product of even powers depends only on the 
exponents and not on their order. 

The proof of this theorem is direct in view of the supposed symmetry of R. 
Now if we impose a requirement of symmetry on a numerical integration formula 
we find that it is possible to compute useful formulas for these regions which 
have not heretofore been available. A numerical integration formula is sym- 
metrical if the set of evaluation points is decomposable into symmetrical sets of 
points, each point in a symmetrical set having the same weight as the others. 


THEOREM 4. Let > a;f(x;) be a symmetrical sum formula approximating the 
integral f f(x)dx for R symmetrical. In order that ¥ a;f(x;) — f f(x)dx = 0 for 
R R 


every polynomial of degree at most 2k + 1 itis necessary and sufficient that the formula 
hold exactly for all monomials f of form 


£,2*1. . .f,, 2m 
where 
ki2>khe--->kn>O0, ki t--- +h Sk. 


We illustrate the application of Theorem 4 by considering symmetrical regions 
in 3-space and indicating a'‘method of obtaining a formula involving at most 34 
points holding exactly for 7th degree polynomials. The exact formula is given for 
a cube. Now in 3-space an arbitrary selection of a point would involve 47 others 
with the same weight. Hence we select points for which the symmetry requirement 
does not yield as many points. In considering the problem we will use (x, y, 2) 
as coordinates. From Theorem 4, we then solve equations only resulting from 
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asking that the symmetrical formula hold exactly for the following seven func- 
tions: 1, x*, x*, x*y?, x®, x+y*, x*y%s? in order to achieve a formula exact for the 
general seven-degree polynomial. For seventh degree in variables with a sym- 
metrical formula over a symmetrical region, the same seven functions suffice. 
Now let the integrals of the seven functions given over R be Jo, I2, Is, I22, Is, 42, 
I222, respectively. We choose points and weights as follows: (x:, 0,0), a: (6); 
(x2, X2, 0), @2 (12); (xs, xs, Xs), 3 (8); (X4, X4, X4), a, (8). The integer in parenthesis 
indicates the number of points in the formula generated by the required sym- 
metry if the indicated letter coordinate is not zero. Since only squares of unknown 
coordinates appear we let u; = x;", u2 = x27, us; = x3? and uw, = x2. 

The equations then are: 


(1) 6a; + 12a. + 8a; +8a,=15 f=1 

(2) 2aiu + 8a2t2 + 8a3u3 + 8a,= 1, f= x 

(3) 2a,;u2? + 8a2u2* + 8a3u;? + 8aueg=—1l, f =x 

(4) 4a2u,? + 8a3u;? + 8aug=I2n f= xy 
(5) 2a,u;* + 8aeus* + 8au;° + Bau =I, f = x® 

(6) 4a.u + 8asu;* + 8aue=Ie f = xy 
(7) Sasus? + Saude = Tenn f = x*y*2*. 


It will be observed that the set constitutes seven equations in eight unknowns 
so that we may hope for one free choice of another equation. To indicate the 
solution of the equation we now form (6)—(7) and (5)—2(6) + (7) to get, 
respectively, . 


(8) 4aous® = Ig. — Io22 
(9) 2ayuy> = Ig — 2D42 + 222. 


Hence if a solution to the equation system exists and the numbers on the right 
of (8) and (9) are not zero then neither u; nor “2 is zere and hence 


(10) 4a2 = (Tae — I222)/u2* 

(11) 2a, = (Ie — 2142 + I222)/t1'. 
Now from (3)—(4) we find 

(12) 2a,u;? + 4aous? = I, — Foo, 
and on substituting from (10) and (11) in (12) we have 


I, — 21 I Ig —TI 
(13) 6 a2 + Lo22 4 i 222 ie eit 
Uy, Ue 









At this point we find it expedient to assume a relationship between u, and u:, 
for example “; = Au. For the cube, later on this shall be taken equal to unity. 
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nc- For a specified 4 we may solve (13) and u; = Aw, for uz and u;, and then by (10) 
the and (11) determine a; and a2. Hence from this point we assume 4, @2, #1, t%2, as 
=m F known, and proceed to find u3, u4, a3, a4. To this end we rewrite equations (1), 
ice. (2), (4), and (7) as 
I 42, 
6); i (1’) 8a3 + 8a, = Io — 601 — 1242 = po 
a, é (2’) 8azuz oa Sag, = IT, ane 2a3;u4; -_ 8a2te = Pi 
wn , (4’) 8a3u;? + 8auu?e = Too — 4a2u?? = Pe 

F (7’) 8azu;* + 8aue = In22 = ps. 


Now if there is a solution to these equations then there exists a quadratic 
function 2? + ¢,z + co such that u; and u, are zeros of the function. Then, how- 
ever, C; and ¢o must satisfy 


Fide ko 


(14) Polo + Pici + P2 = O 
3 Pico + Poti + Ps = O. 


orks 


a 


Assuming such a solution ¢o and ¢; we find u; and u, 


TN 

















} — 1 — Vox? — 45 — 1 + Ver? — 4¢6 
} (15) u3 = . uy = , 
: 2 2 
ye ; and by substitution in (1’) and (2’) we may then find a; and a,. It is not simple 
. a FE to prove that (14) is solvable, nor that u; and u, are positive as determined. It is 
get, quite conceivable that symmetrical regions exist for which some of these wishes 
: are not fulfilled. 
: We now let R be the cube of edge 2, center (0, 0,0) and with edges parallel 
to the axes. Then Io => 8, Te = 8/3, I, = 8/5, Teo = 8/9, I; = 8/7, la = 8/15, 
Ix22 = 8/27 and we find the following solution if we let A = 1, i.e., uw; = us, 
ion 5 _ 1078 ie J! 
1 * 3645 ie 
| O2 = 3645 eee 
: _ 960 — 3V28,798 4, = 200+ 3V28,798 
_— 2726 _ 2726 
4 
5 230 — 774u 
; = Vu, ay: ~~ seein 
Pa re > * 1215 (us — us) 
774u, — 230 
= Vu, Sie. 
RS it + 1215 (us — us) 
d ts, The evaluation points all lie in the cube and the weights are positive, both 


nity. features being desirable for many applications. It is possible to solve the system 
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with other specializations. In case one takes u, = 0 then two 27-point formulas 
evolve, but some of the points of each lie outside the cube as James Clerk- 
Maxwell found in 1877. The 27-point formulas for the sphere will be given in the 
sequel. These too, have points outside the sphere. 

Let us consider an interpretation of the 34-point formula given here for the 
cube. Let f be a function and let the formula be applied to approximate the in- 
tegral over R. If one takes any seventh degree polynomial function agreeing 
with f at the 34 evaluation points, it wili have the integral given by the formula 
no matter how the remaining 120 — 34 = 86 degrees of freedom are used to 
determine the particular 7th degree polynomial since there are 120 terms in the 
general 7th degree polynomial. These 86 degrees of freedom are available to 
achieve better error estimates (but not, of course, a smaller error). The formula 
may be used as a formula for the parallelopipedon by Corollary 1.1. 

Error Analysis. In general, error analysis in higher dimensions is much more 
difficult than for the function of one variable. Nevertheless, it is perhaps of value 
to summarize briefly one type of error analysis. Let 


E(f) = Xadf(e)) - f w(x) f(x)de. 


Let P be a class of functions p such that E(p) = 0. Then P may be considered 
as a linear space since if p, e P and p.¢ P then E(c¢:p: + cop2) = 0. Now for any 
particular f chose pe P and define r by f — p = r. Then 


E(f) = E(o +17) = E(r) = ¥ az(x;) - ff roceyr (era 
R 
This statement holds for all p « P and every appropriate f. 
We want to choose p for a given f so that there is a calculable bound to | E(r)| 
and hence to |E(f)|. First choose p so that p(x;) = f(x;). Then r(x) = 0 and 


E(r) = - J w@rmar so that |E(r)| <b J werax if |r(x) |< b. The associ- 


ated approximation problem is to find pe P so that p(x;) = f(x;) and so that 
max|f — p| is minimized. If this can be done, we have achieved the best error 
bound associated with the last inequality. On the other hand, if one takes pe P 
such that max| f(x) — p(x)| < } for x ranging over R and any ~x outside R and if 


a,>0, La= f wear, 
R 


then . 


|IE(f)\< 2b f w(x)de. 
R 


In applications the class P does not contain all the functions » such that 
E(p) = 0 since these are unknown and essentially unknowable. 

Example. Hammer has devised a code for IBM 650 use which has been checked 
out by Gerald Bartholomew. This code is capable of numerical integration over 
any symmetrical region in E; using 34, 27, 14, 8, or 6 points. For a test calculation 
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the integral of 1/V¥x + 1/¥y + 1/vz, was approximated over the cube with 
center (1/2, 1/2, 1/2) and edge 1. The 34-point formula here described was used 
and the result was 6.04 compared to the exact value 6. 

Conclusion. We have here summarized a few of the useful theorems for generat- 
ing numerical integration formulas from given ones, and illustrated one particular 
approach to obtaining formulas for symmetrical regions. In a sequel we will 
present all the particular formulas we have calculated for symmetrical regions 
including some for circle, square, sphere, cube, hypersphere, and hypercube. The 
method of attack presented for symmetrical regions does not hold much promise 
of effective extension to arbitrarily high-degree polynomials such as Peirce has 
carried out for the circular annulus. However, variations and extensions of these 
methods will be forthcoming. 
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Solution of Poisson’s Equation by Relaxation 
Method—Normal Gradient Specified 
on Curved Boundaries 


Introduction. In the application of the relaxation technique to the solution of 
the Poisson’s equation when the normal gradient of the wanted function, w, is 
specified along a curved boundary, difficulty has been experienced in finding a 
suitable finite difference approximation to the Poisson’s equation at nodes that 
lie adjacent to the curved boundary. It is proposed to present in this article 
results which will be found adequate for tackling any type of curved boundary in 
two dimensions. 

This problem was treated by R. V. Southwell by means of the membrane 
analogy developed by him [1]. It is, however, preferable to base the approximate 
representation of differential equations on finite difference theory, which indicates 
the order of the error involved in the approximation. This has been fully achieved 
in the case of the Poisson’s equation when the boundary condition consists in 





68 SOLUTION OF POISSON’S EQUATION BY RELAXATION METHOD 
(i) the wanted function, w, assuming a given value on a boundary of any 
shape, or 

(ii) the normal gradient of the wanted function assuming a given value on a 
rectangular boundary or a boundary made up of straight lines parallel or per- 
pendicular to the mesh lines. 


The general problem of a curved boundary with the normal gradient of w 
specified thereon, however, has exhibited great technical difficulties in the work 
of L. Fox [2], for example. Nodes in the diagonal directions are introduced into 
the approximate formulae by D. N. de G. Allen [3] in order to tackle this problem. 
In the present treatment approximate formulae are derived by taking into ac- 
count the rate of change of the normal gradient of w along the boundary. 

The formulae developed here lead to simple relaxation patterns. The coeffi- 
cients occurring in these formulae depend on the curvature of the boundary at 
the successive nodes. Their computation at the outset necessarily involves heavy 
arithmetic: The final result can be expected to be correspondingly nearer the 
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3 ° 1 
7 4 8 














Fic. 1. Neighbouring nodes numbered with reference to Node ‘0’. 


correct solution. Close adherence to the specified boundary condition is ensured 
as both the curvature of the boundary and the rate of change of the normal 
gradient of w along the boundary are taken into account by the formulae de- 
veloped here. 


Symbols and Notation 


w the function which is to be evaluated in a given region 

W a function given in the region of integration 

k a function given on the boundary of the region of integration 

v? the Laplacian operator 

0 the node at which a residual formula is required 

Fy residual at 0 

1,2,---,8 Neighbouring nodes that surround 0; these are to be identified 
according to the scheme shown in Fig. 1 

x9 Cartesian co-ordinates with respect to suitable axes of reference. 

The position of the origin and the orientation of the axes are 

different in different contexts, and are defined where required 
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mesh size 

radius of curvature with the proper sign, plus or minus 
distance measured along the boundary 

an angle defined as in Fig. 4 or Fig. 5 

a ratio defined as in Fig. 4 or Fig. 5 

tan a. 


The Problem. Determine the value of a function w inside a region R given that 


(i) V’w = W a given function, at points inside R, 
(2) dw/dv = k on the boundary of R. 


Here 0/dv denotes the rate of change along the outward-drawn normal to the 
boundary, and & is a function given along the boundary. 


E 
Edge NodeA 
Fic. 2 Fic. 3 


Network and Types of Nodes. For solving the above problem, the region will 
be covered by a network of square meshes of suitable mesh size as shown in 
Fig. 2, and the value of w will be ascertained at the nodes of this network by the 
method of relaxation. The diagram of the network is shown separately in Fig. 3, 
as this diagram will engage our attention throughout the process of relaxation. 
On examining this diagram, we find that there are three types of nodes in the 
network. 


(a) For most of the nodes (i.e., nodes not belonging to the two types discussed 
below) there are four neighbouring nodes, which can be numbered 1, 2, 3, and 4 
in accordance with the scheme set out in Fig. 1. These nodes will be called 
“interior nodes.” 

(b) For nodes which occur on the edge of the network, one of the neighbouring 
nodes 1, 2, 3, or 4 is missing; e.g., node E for which there is no neighbour which 
can be marked with the number 4. Such nodes will be called ‘‘edge nodes."’ 
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(c) For nodes which occur at a corner of the network two of the neighbouring 
nodes are missing; e.g., node C for which there are no neighbours which can be 
marked 2 or 3. These nodes will be called ‘‘corner nodes.”’ 


In the case of an interior node it is easily derived [3] that the finite difference 
approximation to the Poisson’s equation is 


(3) W2Wo = wi + we + ws + ws — 4w0 + O(h'), 


where the suffixes denote the nodes at which the values of the concerned functions 
are to be taken. This result is obtained on assuming that the wanted function 
possesses a Taylor’s series at the point O, with a circle of convergence large enough 
to enclose nodes 1, 2, 3, and 4. 

The above approximation is of no avail at the edge and corner nodes, as at 
these nodes one or two of the neighbouring nodes, 1, 2, 3, and 4, are missing. The 
boundary condition which w has to satisfy, has to be utilized in order to eliminate 
the values of w at the missing nodes from equation (3). This is achieved as follows. 

Edge node. Taking an edge node 0 let it be assumed that node 1 is missing 
(Fig. 4). Let it be assumed, again, that the wanted function, w, possesses a 





—— kr 





. 


Fic. 4 


Taylor’s series at 0, with a circle of convergence large enough to enclose nodes 2, 
3, and 4. 
Let the series be 


(4) w= wot px + qy + bx + 2fxy + cy 
+ ayx* + 3bix*y + cry? + diy? + ---, 


where x, y are cartesian co-ordinates with origin at 0 and x-axis along the straight 
line 3-0-1. 
Since V-w = W at any point in R, 


(5) 2(b + c) = Wo. 


Let the straight line 3-0-1 cut the boundary at A. The distance OA is less than 
or equal to h. Let OA = th, where 0 < — < 1 (vide Fig. 4). Let the normal at A 
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to the boundary make an angle a with the x-axis (i.e., 3-0-1), in the anticlockwise 
direction. 
The given boundary condition states that 


(6) (@w/dv), = ka. 
Since dw/dv = k all along the boundary, it can be derived that 
(7) (0*w/dsdv)4 = (dk/ds)a, 


on assuming that k possesses a derivative along the boundary at A. [0/ds denotes 
the rate of change along the boundary. } 

It is easy to express (@w/dv),4 and (0*w/dsdv), in terms of the coefficients in 
the expansion (4) since along the curve, 


(8) dw/dv = dw/dx cos (x, v) + dw/dy cos (y, v) 
and 


(9) &w/dsdv = d°w/dxdy[cos*(x, ») — cos*(y, v) ] 
+ [6*w/dy* — d*w/dx*] cos (x, v) cos (y, v) 
— 1/p[dw/dx cos (y, ¥) — dw/dy cos (x, v) ]. 


Here p is the radius of curvature given by ds/dy in the usual notation. The 
symbols (x, v) and (y, v) indicate the angles between the x- and the »-directions, 
and the y- and the »-directions, respectively. 

By differentiating the series (4) successively and substituting x = th, y = 0 
in the results, the following approximate values are obtained at the point A. 


(10) (dw/dx), = p + 2th + O(h*), 
(11) (dw/dy)a = g + 2fth + OCF’), 
(12) (d*w/dx*), = 2b + O(h), 

(13) (@’w/dxdy), = 2f + O(h), 

and 

(14) (0°w/dy*)a = 2c + O(h). 


Combining (6), (8), (10), and (11) we obtain 
[A] ka = pcoosa + gsina + 2th(b cosa + f sina) + O(h’). 
Similarly combining (7), (9), and (10) through (14) we obtain 


[B] (#) mp 4 gO + 2 f c08 2a + 6 — b sin acosa} + O(h). 
ds/A Pa Pa 
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On multiplying [A] by seca and [B] by £¢h tanasec 2a and taking the 
difference, 
(15) kaseca — th tan asec 2a(dk/ds) 4 
= p(1 + éhsin a tan a@ sec 2a/pa) + g(tan a — th sin a sec 2a/pa) 
+ 2bth — (c — b)th tan a tan 2a + O(h). 
On substituting the co-ordinates (0, h), (—h,0) and (0, —h) respectively in 


the series (4), the values of we, w3, and w, are obtained. From these it can be 
deduced that 


(16) We + ws — 2wo = 2ch® + O(h') 


We—- Wy = 2hq + O(h') 
and 
W3 — Wo = — ph + dh? + O(h*). 
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Fic. 5 


Combining (5), (15), and (16) to eliminate the unknown coefficients, we obtain 
an elegant result, 


1 + m? 
1 — m? 





dk 
(17) Woh? (1 + 2é — - 2hkaV1 + m? + m + 2(¢ ) th?m 
m* ds/JA 





1 h mN1 + m? 
= wi (1 — m + 267+ = amin) 


h m?N1 =e 
+ imi i + ¢-— 
min 


1+ h mite) 





+w (usnanits iy ie ye 


Fad haem + at 
pA i= 





— 4wo (1 + &5 ) + O(f'), 


where m = tan a. 





the 


(h*). 


y in 
1 be 


fain 








SOLUTION OF POISSON’S EQUATION BY RELAXATION METHOD 73 


The first member of this equation can be calculated at every one of the edge 
nodes, as soon as a suitable square net is fitted to the region of integration, since 
this expression depends only on the given functions W and k, and the character- 
istics of the boundary curve. Similarly all the coefficients of the w on the right 
hand side can be determined at the outset. 

Corner node. The method of obtaining a residual formula at a corner node 
involves a few geometrical constructions, which are indicated in Fig. 5 for the 
case in which nodes 1 and 2 are missing. Taking our stand for a moment at 0’ in 
Fig. 5, we can derive as we have done above, that 


(18) Wee [1 +26 —|- V2hkaV1 + m+(4 EY em * ha 








1 — m? 
+ m? h at) 
= ws (1 — m+ 264 i-w “% i-s 
42 (14 h mvNi+m m) 
W7 a i- 





eat (1 +m +t - h wt) 
m? V2p, 1+m 


1 h Vi + m? 
— tue (1 +6 tem te 
ee. 1—m 


The positions of 0’ and A are to be noted carefully. Again 





=) + OCF). 


0-3 = 0'-7 = 0-4 = h/V2. 


Hence (18) is obtained from (17) by merely replacing h by h/ V2 and giving suffixes 
to w in accordance with Fig. 5. It is to be noted that ¢ > 1 here. 
Because of Poisson’s equation, 


h2 
(19) — 9 We = Wo + Ws + Wr + ws — 40 + OC). 
On eliminating wo between (18) and (19) we obtain 
ry Oe Vanes Vid (*) tn Lt 
(20) Wo Zé 2hkaVi + m* + ds Patek WE 
‘<< th ee) 
s(- m + fF V2. 1 — m2? 
(1 1 + m? th = ee 
=m << V2p4 1 — m?* 














+ wa(m + ¢ - 


1+m th =3 2?) 
—m 2, 1—m* 








1 + m? th mette) 
(1 +€é ee + Vans "ae + O(h'). 
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The important points about this formula are: 


(i) The value of W to be taken is the value at 0’ and not that at 0. 
(ii) A change in the assumed value of w at node 7 will affect the residual at 0 
but not vice versa. 


Illustrative Example. The above formulae were applied to calculate the values 
of w = logr = } log (x? + y’) within the region BCDEFB shown in Fig. 6, 
assuming dw/dv to be given all along the boundary. It is well known that log r 
satisfies Laplace’s equation in any region of the plane, provided the origin is 
outside the region. Hence W = 0 in equation (1). 

The curved boundary FB is part of the ellipse 


x2/4 + 2 = 1. 


BC and EF lie along the x- and y-axes respectively, while CD lies along the line 
x = 3, and DE along the line y = 2. 

The network adopted is rather coarse with 4 = 0.25 only. There are eight 
nodes besides F and B, which lie adjacent to the elliptic boundary. They are 
numbered in order in Fig. 6. Of these nodes, NV; and N; are corner nodes and the 
rest, edge nodes. Hence equations (17) and (20) are directly applied to these 
nodes. Table 1 summarizes all the data required at these nodes for calculating 
the coefficients in equations (17) or (20), and also gives the corresponding residual 
formulae. 


TABLE 2 
Nodes Residual 
F Fo = W, oe We2— 0.250 
B Fo = Wi +.wWe — 0.125 
Cc Fo = We + wz + 0.083 
D Fo = wz + ws + 0.096 
E Fo = Wt w, + 0.125 


At nodes F and B the missing nodes are eliminated from equation (3) by 
making use of the values dw/@x and dw/dy which are both known there. The same 
is the case at nodes C, D, and E. The residual formulae for these nodes are given 
in Table 2. 

At all intermediate nodes on BC, CD, DE, and EF, one of the surrounding 
nodes is missing and this is eliminated from (3) as either dw/dy or dw/dx is known 
at these nodes. (Alternatively the same result can be obtained from equation 
(17) by setting § = m = 1/p = 0.) The residual formula at a node on BC, for 
example, will take the form 


Fo = Wi + 2wWe + W3 — 4wo + 2hko = O(h'). 


It may be noted that the problem as set up so far is indeterminate in the sense 
that any constant value can be added to a solution that satisfies the conditions 
of the problem. For 


ll 


if Vw 
and dw/dv = 


W inside a region R, 


k on the boundary of R, 
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it is clear that 
V?(w +c) = W inside R, 





























































































































and 
at 0 0/dv(w + c) = k on the boundary, 
where c has the same value at all points. This property is reproduced in the set 
lues of residual formulae derived above. Every one of these formulae has the feature 
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Fic. 7. The End: Values of w multiplied by 10* are shown. Final residuals are given 











78 SOLUTION OF POISSON’S EQUATION BY RELAXATION METHOD 


that the sum of the coefficients of the unknown w is zero. Hence if a set of values 
of w give rise to a certain set of residuals, then the set of values (w + c) where c 
is the same for all nodes, will also give rise to the same set of residuals. It is, 
therefore, necessary to give the value of w at some point in order to make the 
solution unique. We make use of the fact that log r = 0 at F for this purpose. 

Starting with the initial guess w = 0 at all nodes, the solution exhibited in 
Fig. 7 was obtained by the method of relaxation. The initial residuals are shown 
in Fig. 6A. During the process of relaxation of the residuals, negative values were 
obtained for w at nodes like F, N,, etc. (vide, Fig. 6). When all the residuals had 
become negligible, a suitable value was added at all the nodes so as to bring the 
value of w at F to zero. 

The values of log 7 for select values of r as taken from Fig. 7 are compared in 
Table 3 with the true values taken from a standard table of Mathematical 








TABLE 3 
Log r Log r 

by re- true by re- true 

r laxation value Node r laxation value 
1.00 0.000 0.000 Ni 1.0308 0.030 0.030 
1.25 0.221 0.223 N2 1.1180 0.110 0.111 
1.50 0.402 0.405 N3 1.2500 0.221 0.223 
1.75 0.556 0.560 N, 1.4142 0.345 0.347 
2.00 0.689 0.693 Ns 1.4577 0.373 0.377 
2.25 0.806 0.811 Nes 1.6771 0.513 0.517 
2.50 0.911 0.916 N; 1.8200 0.594 0.599 
2.75 1.006 1.012 Ns 2.0156 0.696 0.701 
3.00 1.093 1.099 D ¥13.0000 1.279 1.282 


Functions, see F. Castle [4]. The agreement is remarkable in spite of the network 
being coarse. 

The author is indebted to D. N. de G. Allen of the Imperial College of Science, 
London, for the interest shown by him in the publication of this article, when he 
visited the University of Roorkee, India. 

Messrs. C. M. Ganesan and G. Natarajan have been of valuable help to the 
author in the preparation of the typescript. 

R. V. VISWANATHAN 
Central Water and Power Commission 


Government of India 
New Delhi, India 


- 1. = SOUTHWELL, Relaxation Methods in Theoretical Physics, Oxford, at the Clarendon 
ress, 1946. 

2. L. Fox, ‘“‘Numerical solution of elliptic differential equations when the boundary conditions 
involve a derivative,’’ Roy. Soc., Phil. Trans., London, v. 242, 1950, p. 345-378. 

3. D. N. pE G. ALLEN, Relaxation Methods, McGraw Hill Book Co., Inc., New York, 1954. 

4. F. Caste, Five-Figure Logarithmic and Other Tables, Macmillan and Co., London, 1910. 
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TECHNICAL NOTES AND SHORT PAPERS 
Table Pertaining to Solutions of a Hill Equation 


1. Introduction. The electronic digital computer of the Graduate College of 
the University of Illinois has been employed to determine the characteristic 
exponent and a quantity related to the phase and amplitude functions for solu- 
tions of a Hill equation of the form 


@Y/d? + (A + Boos 2t + Ccos 4t + Dcos 6f)Y = 0. 
Results were obtained for the following values of the coefficients: 


D=0 

C: —0.5(0.1)0.5 

B: 0(0.2)0.8(0.1)1.6(0.2)2.0 

A: -—0.5(0.1)—0.2(0.05) —0.15(0.03)0.15 (0.05)0.2(0.1)0.5 ; 
and 

D: —0.05(0.05)0.10 


C=0 
B: as before 
A: as before. 


2. The Quantities Determined. The quantities which were determined can be 
defined by reference to the matrix M(r) which relates the value and slope of a 
general solution at ¢ = + + = to the value and slope at ¢ = r: 


steed min (or) (= mba } 
= T = ° 
Y'(r + ™) Y'(r) Mx M2 Y’(r) 
(i) The first quantity of interest is the invariant 
4 Tr M = cosh typ, 
where yu is the characteristic exponent of the Floquet solutions [1 ]. In cases such 
that — 1 <}TrM <1, corresponding to stable solutions of the differential 


equation, it is convenient to write 


4 Tr M = cosc. 


In terms of the phase-amplitude representation of such solutions, Y(t) = w(t)e+** 
= w(tetlel™++¥(] (with w(t) and W(t) periodic in t, mod z, and with w’ a 
constant), o is interpretable as the advance of the phase when ¢ increases by r. 

(ii) The second quantity of interest is the matrix element M,2(r), which may 
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be identified with a . In terms of M,2(r) the complete matrix may be written 
T 
cos ¢ — 4M‘. Mi. 
Mr) = 1(M",.)? +n? : 
(r) a i( 12) + sin" ¢ cos ¢ + 3M’ 12 
My 


while the amplitude function w(r) may be taken as proportional to VMi:. Because 
of the properties mentioned it was felt that, in the interests of brevity, it would 
suffice to characterize the solutions by tabulating Mi2(r). Since this quantity is 
an even periodic function, it is only necessary to tabulate values in the range 
0<¢7€ 2/2. 

3. Method of Computation. For each set of coefficients A, B, C, D two solu- 
tions of the differential equation were obtained in the course of the computation. 
The solutions were calculated by direct integration of the equation, using the 
Runge-Kutta procedure with a step size h = 2/64. The initial conditions imposed 
on these solutions were 

Y,(0) =1, Yi’) =90, 


and 
Y2(0) = 0, Y,’ (0) = 1. 


Because of the symmetry of the problem, the integration was restricted to the 
range 0 < t < 2/2. 
The quantities which were tabulated were then determined by the equations 


(1) 1/2 Tr M = Y,(x/2) Vo’ (x/2) + Vi! (4/2) V2(x/2), 
(2) Mi2(r) = 2{ Vo(m/2) Vo! (e/2)CVil(r) P — Vilw/2) Vi’ (w/2)LY2(r) FP}. 


(As is shown in the introduction to the completed table, these relations may be 
readily established by an appeal to symmetry which permits matrices which 
carry solutions from 2/2 to z and from ¢ to 2/2 to be Written in terms of those 
which apply to the range 0 to z/2 and 0 to ¢.) The matrix element M2 was thus 
computed for the nine values of r: 0(7/16)x/2. In those cases for which 


1 
—1<3TrM <1 the quantity © = cos! (3 Tr M) was also computed, the 
T T 


inverse cosine being taken to lie in either the first or second quadrant. 

4. Estimated Accuracy. The chief error in the cgmputation was regarded as 
arising from truncation error in the integration. As is outlined in the introduction 
of the table, this error was estimated (i) by considering the absolute error in a 
single integration step to be h' YY/120, and (ii) by repeating representative calcu- 
lations with a step size h’ = 4h. These two methods of estimating the error 
led to the conclusion that the tabulated values of cos ¢ and of My: are in error 
by no more than 4 X 10~-*, or by several units in the last decimal place retained. 
In the case that |cos ¢| is very close to unity, a small error in cos o will of course 
be magnified in the evaluation of o/z. 
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5. Availability of Tables. The table described in the present article has been 
deposited in the file of Unpublished Mathematical Tables which MTAC main- 
tains. This material may be made available on loan to interested people [2]. 


G. BELFORD 
University of Illinois 


Urbana, Illinois 
and 
Midwestern Universities Research Association 
L. Jackson LASLETT 
On leave of absence from 
Iowa State College 
Ames, lowa 

J. N. SNYDER 
University of Illinois 
Urbana, Illinois 
and 
Midwestern Universities Research Association 


The authors were assisted by the National Science Foundation and the Office of Naval Research. 


yy A | eee & G. N. Watson, A Course of Modern Analysis, Cambridge Univ. Press, 
1946, p. 412 ff. 

2. G. BELForp, L. Jackson Lastett, & J. N. SNyDER, Table Pertaining to Solutions of a Hill 
Equation [MTAC, this issue, Review 59, p. 110]. 





Rational Formulae for the Production of a Spherically 
Symmetric Probability Distribution 


Von Neumann [1] has described the following method for producing a uniform 
distribution on the circumference of the unit circle: Select two random numbers 
*1, X2 from a uniform distribution on (—1, +1). Reject the pair if x;? + x? > 1, 
and re-select until a pair is obtained satisfying x; + x.* < 1. Then, by the trigo- 
nometric double-angle identities, 


x? - x? 2x1X2 
Be ’ y SS ee 
x + x? xe + x? 


are the co-ordinates of a point having the desired distribution. 

The following is an extension of this method for points distributed on a sphere 
uniformly with respect to area. Select four random numbers xo, x1, x2, X3 from a uni- 
form distribution on (—1, +1). Reject the quadruplet if x? + x:* + x2 + x; > 1, 
and re-select until a quadruplet is obtained satisfying xo? + x;? + x2? + x;* < 1. 
Then 


_ __ 2 (%1%3 + xore) _ _2(%a%3 — xox) ga et ae — ait — 27? 
xo? + x? + xo? + x3?’ / xo? + xy? + x2 + x3?’ xe + x + x2? + xi? 











are the co-ordinates of a point uniformly distributed over the surface of the unit 
sphere. 
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A proof depends on the fact that a non-zero quaternion q = xo + xii + xoj + x3k 
effects the rotation v — qvq-' of vectors v, and that conversely every rotation 
can be so produced by some quaternion. 

Select q from a uniform distribution on the set |q| < 1 as above. Then qkq~ 
will have on the surface of the unit sphere a distribution invariant under rotations 
because, if p is any other non-zero quaternion, then 


p(qkq~*)p™' = ( in] a) k (2 a) ; 


but multiplication of q by a quaternion of unit norm is a rotation of quaternion 
space and hence preserves the distribution from which we selected q. 

But there exists on the surface of a sphere only one distribution invariant 
under all rotations, namely, the uniform distribution. 

Now let qkq™ = xi + yj + 2k and we have the above formulae. 

The method extends to » dimensions and to more general spaces acted upon 
by a compact group of transformations. 

J. M. Coox 


Argonne National Laboratory 
Lemont, Illinois 


This. work was performed under the auspices of the U. S. Atomic Energy Commission. 


1. JOHN VON NEUMANN, “Various techniques used in connection with random digits,’ NBS 
Applied Mathematics Series, No. 12, U. S. Gov. Printing Office, Washington, D. C., 1951, p. 36-38. 


On the Numerical Solution of Integral Equations 


For the numerical solution of linear integral equations we must of course 
replace all our variables, both dependent and independent, by discrete variables. 
In choosing the grid for the independent variables we should naturally be influ- 
enced by the storage facilities of our calculating machines. One purpose of this 
note is to emphasize the obvious fact that we should not be too strongly influenced 
by the storage facilities since we may then take longer than necessary for the job. 
The following technique is formulated in terms of finding the inverse of a kernel, 
K(x, y), but a similar technique could be applied for solving a single integral 
equation. The idea seems too simple to be new but I have not met it before, and 
it may have some applications. 

First choose a wide grid. This will convert our kernel into a small matrix, M,. 
Invert it. Next choose a less wide grid, and hence a larger matrix, M2. Use the 
result of stage 1 to construct an approximate inverse of M; (by smoothing the 
elements of the inverse of M; so as to obtain a matrix of the right size). This step 
will usually be justifiable if K(x, y) is a smooth function. Starting with our ap- 
proximate inverse of M, we may obtain a more exact inverse by means of a well 
known iterative procedure known to be rapidly convergent when we start with 
a good approximation (I. J. Good [1] and H. Hotelling [21]). We can next choose 
a still narrower grid and proceed as before. The whole process can be stopped 
when adjacent entries in our matrix are numerically close. 
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In one-dimensional problems of mathematical physics the word “adjacent” 
can here be usually interpreted positionally on the matrix itself, if the rows and 
columns are in an appropriate order. Otherwise the word shoulc be interpreted 
in terms of the physical application. For example, in the inverse of the matrix 
L, treated by W. L. Wilson [3], which arises from a two-dimensional physical 
problem, the elements at the positions (4,14) and (7,19) may be regarded as 
adjacent because they correspond to pairs of points that are close together in the 
physical problem. This question of the definition of “‘adjacent”’ is related to that 
of how to extrapolate from a matrix to a larger one. It seems difficult to formulate 
any general rule for extrapolation. For the present each case would have to be 
treated on its own merits by the exercise of judgment, as in the reference 
cited. 

I. J. Goop 
25 Scott House 


Princess Elizabeth Way 
Cheltenham, England 


a7 J. Goon, “Bounded integral transforms,” Quart. Jn. of Math., Oxford, v. 1, 1950, p. 
185-190. 


2. HAROLD HOTELLING, ‘‘Some new methods in matrix calculation,"’ Annals Math. Stat., v. 14, 
1943, p. 1-34. 
3 


. W. L. Witson, Jr., Tables of Inverses to Laplacian Operators over Triangular Grids, UMT 
Fite [MTAC, this issue, Rev. 53, p. 108]. 


A Rotation Method for Computing Canonical Correlations 


Given two sets of variates, x;(i = 1,2 --- p), y;(j = 1,2 --- g) with p 2 gq, 
Hotelling [1] has shown that it is possible to find linear transforms u,, v; of the 
«’s and y’s respectively with the properties 

1. var (u;) = var (v;) = 1 
2. cov (u;,um) = 0, i#k 
cov (v;,m.) = 0, 7 #1 
3. cov (ui,0j)) = 0, 4 j. 
The variates u and wv are called canonical variates, and the correlations 


pi(t = 1,2 --- g) between corresponding variates u;, v; are called canonical 
correlations. 


The same problem may be stated in terms of matrix algebra. Suppose that 
the dispersion matrix of the x’s and y’s, considered as a single vector variate, is 


A C 

Cc’ B 
where A has p rows and columns, B has g rows and g columns, and C has p 
rows and g columns. The whole matrix is symmetric and positive definite. We 
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seek matrices U of » rows and columns and V of g rows and columns such that 
UAU’ =I VBV’'=I 
UCV’=R 
where R = (rj), Tis = pi, Ti3 = 0,4 F 7. 
Hotelling showed that the p,? were the roots of the equation 
|C’ AC — »B| = 0 


and that the coefficients of v; could be found as the solutions of the consistent 
homogeneous equations 


(C’A“C — p;B)y = 0; 


the coefficients of u; are then given by A C y. This procedure clearly presents a 
formidable computational problem if p and g are at all large. If a high-speed 
computer is available, the problem is greatly reduced, since the various stages 
form standard computational problems; A~! C and C’ A-' C can be calculated by 
a process of pivotal condensation as described by A. C. Aitken [2], and the roots 
and vectors derived from (C’ A~! C — \B) can be readily obtained from the latent 
roots and vectors of B-'C’ A-'C, or from those of the symmetric matrix 
K-' C’ A-'C K’-", where K is a triangular matrix such that B = K K’. Never- 
theless, several fairly complicated stages of computation are involved, and their 
linking together is not altogether straightforward ; in addition, a fair amount of 
intermediate storage capacity is called for. The purpose of this note is to describe 
a more direct method of computation which avoids some of these difficulties. 
The method suggested here is similar to a well known rotation method for 
finding the latent roots and vectors of as ymmetric matrix (see Householder [3 ]). 
We wish to find matrices U of » rows and columns and V of g rows and columns 


such that 
U Oo A C ee). 4. B 
0 V . = 0 vJ \R I 


Triangular matrices H and K can be found such that A-= H H’ and B = KK’ 


(computationally, this is a straightforward and accurate process); writing 
U H = S and VK =T, we obtain 


S 0\/H" 0 \/A C\/(H' 0 yeh f.2 
o T/\o K-*J\c B/\o x=/)\o rr)” \R' 1 
S 0 I H“CK’\ (s’ 0) _ (IR). 

o tT) \k-c’H I i a 


In this form it is seen that we require two orthogonal matrices S and T such that 
they together transform the matrix H-! C K’— of » rows and g columns to the 
diagonal form R. We can now show that S and T can be found iteratively as the 
products of sequences of elementary rotation matrices. 

Writing for brevity H-' C K’' = D, consider the matrix 


D,; = S; D =, 


or 
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where 

cos6; sin@, 0 te 

—sin@; cos@é,; 0---0 

Ss; = 0 0 1---0 

0 0 0---1) 

cos 6, sin 02 Ber 8) 

—sin 6, cos@, 0---0 

T= 0 0 1---0 

0 0 1} 





It is found that 


10 yC2+-d2151C2+-d1261S2+2251S2 © —b4161S2— de1S1S2+di26162 +225 162 dCi +-dogs1. . . 
D,= — dy381C2+-de1C162— d1251S2+d2201S2 G181S2— d21C1S2—y281C2+-de2eic2 = — 1301 +-des5).. . 
d3,C2+ds252 —dy3iS2+dy2C2 dz; 
etc., 
where 51, C1, S2, Cz are written for sin 6;, cos 6, sin 42, cos 62, only elements in the 
first two rows and columns being affected. 
We now choose 6; and @2 to satisfy the equations 


— dyyC1S2 — deiSyS2 + dy2CyC2 + do2S1c2 = 0 
— dy81C2 + dei€iC2 — dy251S2 + doeeis2 = 0. 


With these values of 6; and 62, the transformation causes two of the elements of 
D, to vanish. 

A similar pair of transformations can now be applied to D, to remove two 
other elements. This will in general restore to non-zero values the elements just 
dealt with; but it will be seen that each pair of transformations dealing with 
elements d;; and d;; leaves unaltered the sum of squares of the remaining non- 
diagonal elements. It can be shown that, on successive iteration, the sum of 
squares of all the non-diagonal elements will tend to zero and the matrix itself to 
diagonal form. The required orthogonal matrices are the products of all the suc- 
cessive rotation matrices. 

The above argument holds good when p = q, so that D is a square matrix. 
If p > q, there will be x variates that have no corresponding y variates. These 
can be dealt with by a transformation of the type D, = S; D which annihilates 
a single element and retains the sum of squares of the remaining non-diagonal 
elements as before. 

The angles 6; and @, are obtained from the equations given above. These are 
found to give 

2 (dirdax + do2d12) 

a? 1, oi d* >, + d* 2 — G22 


2(dirdi2 + dood21) 
iy, + Pa — Pye — oe 


from which the sines and cosines can be derived without difficulty. 


tan 26; = 





tan 20, = 
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A point of practical importance is that the d’s (which are, in fact, correlation 
coefficients) do not exceed 1 in magnitude; the same is true of the r’s. For this 
reason, the iterative scheme can easily be programmed using fixed-point arith- 
metic, though double-precision working is frequently needed in the calculation of 
the sines and cosines. 


M. J. R. HEALY 
Rothamsted Experimental Station 


Harpenden 
Herts., England 


1. H. HOTELLING, ‘Relations between two sets of variates,’’ Biometrika, v. 28, 1936, p. 321-377; 

2. A. C. AtTKEN, “The oy with — of a certain triple product matrix,’ 
Roy. Soc. Edinburgh, Proc., v. 57, 6, p. 172-18 

3. A. S. HOUSEHOLDER, Principles of Numerical Analysis (Sect. 4.115), McGraw-Hill Book 
Co., Inc., New York, 1953. 


Polynomial Approximations to Bessel Functions of Order Zero 
and One and to Related Functions 


In the course of preparing a programme for a large digital computer, the 
following formulas, which are of use in the calculation of Bessel functions on such 
machines, have been obtained. The methods of obtaining the approximations are 
those described by Lanczos in NBS AMS9 [1], with trivial modifications (a 
complete set of references is given in [1 ]). The results enable Jo(x), Ji(x), Yo(x), 
xY,(x), and Jo(x) to be calculated for all values of x and Ko(x), Ki(x), and 
R(x) to be calculated for values of x exceeding unity; where we define 


Jo(x) = f * Jo(t)dt 


and 
Ry(x) = f * Ko(t)dt. 


We define auxiliary functions such that 


Jn(x) = V2: (P.0@ cos (« - = - ) — Qn (x) sin (5 = = “ *)) 


n= 0, 1, 
Y,(x) = 2. (>. (x) sin (.-% -->- *) + Qn(x) cos (« _ > _ *)) 

n= Q, 1, 
9.(x) = (4.0) log = ae Ya(s)) n = 0,1, 
K,(x) = oF Gale), n= 0,1, 


Jo(x) =1—- vi (Po) cos (« + *) — Qo(x) sin (= + *)), 
Rix) = 5 - oy Z a, 
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ion Then with the maximum error stated in brackets in each case, and provided 
his 0Ozc<t<gi 

th- 

of Jo(4t) = 1.00000,00000 — 3.99999,987217 + 3.99999,73021# 


— 1.77775,60599t® + 0.44435 ,84263t — 0.07092,534922" 
+ 0.00767 ,71853#2 — 0.00050,144152" (10 x 10-”) 


1 
7 Js) = 1.99999,99998 — 3.99999,99710f + 2.66666,605441 


- — 0.88888,396492° + 0.17775,82922t — 0.02366,167731 


i + 0.00220,69155¢ — 0.00012,89769 (3 x 10-#) 

1. (4 
+ =P (*) = 0.39894,22793 — 0.00175,30620F + 0.00017,34300# 

2x —\t/— — 0,00004,87613¢* + 0.00001,73565¢* — 0.00000,37043" 

(12 x 10-™) 
1 lta, (*) = — 0.01246,69441 + 0.00045,64324" — 0.00008,69791# 

es tN2r“"\t/ + 0,00003,42468/* — 0.00001 ,420781* + 0.00000,323120" 
- (25 x 10-¥) 
are i 4 
(a Pr (*) = 0.39894,22819 + 0.00292,18256# — 0.00022,32030¢ 
- 2x —\#/ + 0.00005,80759#* — 0.00002,00920¢* + 0.00000,424149" 
a (15 x 107) 


4 
1 ilo, (*) 0.03740,08364 — 0.00063,90400# + 0.00010,64741# 
tV2x \t/ — _ 9,00003,987082* + 0.00001,62200¢* — 0.00000,36594" 
(25 X 10-) 
Go(4t) = — 0.57721,56649 — 1.69113,74142f + 3.69113,887931 
— 2.23311,02234° + 0.66943,21484/* — 0.12141,875612 


+ 0.01489,99271# — 0.00135,08487#* + 0.00008,91322#'* 
(6 X 10-") 


4tg,(4¢) = 1.00000,00004 — 0.61772,53972f — 10.76454,727244 
+ 11.62078,91416# — 4.91052,91148/* + 1.14180,330127° 


A, — 0.16910,81720#% + 0.01699,218762" — 0.00102,66368'* 
(6 X 10-”) 
= mo Ge (*) = 1.25331,41373 — 0.15666,41816¢ + 0.08811,127827 
: t — 0.09139,09546 + 0.13445,96228 — 0.22998,503282° 
» i, + 0.37924,09730# — 0.52472,773312? + 0.55753,683672* 


— 0.42626,32912# + 0.21845,18096¢ — 0.06680,9767274 
+ 0.00918,938302" (10 x 10-) 


Fo G, (*) = 1,25331,41373 + 0.46999,27013¢ — 0.14685,82957# 

2 \t/ 4 0.12804,26636# — 0.17364,31637H + 0.28476,181499 
— 0.45943,421172* + 0.62833,806817 — 0.66322,95430# 
: + 0.50502,38576 — 0.25813,03765° + 0.07880,00118¢" 
| — 0.01082,41775# (10 x 10-") 
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“2 0(8t) = 7.99999,99990 — 42.66666,64204/? + 102.39999,00866/4 

. — 130.03159,01993¢® + 101.13454,38222¢8 — 52.95243,717452° 
+ 19.89838,57672¢? — 5.60245,72363¢ + 1.20146,774492'° 
— 0.18673,71001#* + 0.01624,754657/ (10 x 10-) 


mE P, (®) = 0.79788,45600 — 0.01256,42405# + 0.00178,70944 

© t — 0.00067,401487° + 0.00041,00676r2 — 0.00025,439552"° 

+ 0.00011,07299#2 — 0.00002,26238# (8 x 10-9) 

1 m (°) 
t °’%r Q t 
7 1 
+0.(;) 
V an 


— 0.06233,47304 + 0.00404,03539# — 0.00100,89872# 
+ 0.00053,661692* — 0.00039,9282528 + 0.00027,550372" 
— 0.00012,70039¢" + 0.00002 ,684827* (8 x 10-) 


1.25331,39163 — 0.78323,44963t + 1.25733,120337 

— 3.09054,43850# + 9.02560,45356# — 25.43912,19592¢° 

+ 60.46288,82856¢* — 112.80726,52384#? + 158.83274,70627#° 
— 163.74821,02377 + 119.22659,27008/° — 57.88900,965157" 
+ 16.78876,587872 — 2.19748,244497"* (25 x 10-” in Ko). 


An attempt to derive similar formulas for Io(x), I1(x), and I(x) by the same 
methods proved unsuccessful, since these functions cannot be defined by a 
differential equation and a boundary condition at infinity. 


A. J. M. Hitrcucock 


United Kingdom Atomic Energy Authority 
Industrial Group Headquarters 

Risley, Warrington, 

Lancashire, England 


1. NBS Applied Mathematics Series, No. 9, Tables of Chebyshev Polynomials S,(x) and C,(x), 
U. S. Govt. Printing Office, Washington, D. C., 1952, Introduction. 


An Iterative Method for the Solution of Linear Equations 
Based on the Power Method for Proper Vectors 


Introduction. A computing machine program for obtaining the largest proper 
value and proper vector of an (m + 1,” + 1) matrix A by the power method 
(i.e., by the iteration X,,;; = AX,) may be used for solving a system of m linear 
equations in m unknowns. In the following note a metfiod of setting up the itera- 
tion together with a simple criterion for convergence of the iteration are given. 

Method. Let B be an (m + 1, + 1) matrix with real elements and real proper 
value k such that || is greater than the moduli of all the remaining proper values 
of B. If X is a non-exceptional non-zero column vector, X = col (£1, £2, - - +, En, n41); 
with real components, it is well known that the iteration Xo = X, X; = BXo, 

-+, Xv41 = BX, is such that X,,, = RX,. (There are exceptional choices of X 
for which the iteration will not converge to the largest proper value and vector; 
but round-off error during the iteration will usually make the X, non-exceptional 
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and convergence will actually take place, though possibly quite slowly.) To keep 
the components within size it is usual in a practical computation to modify X, 
after each iteration by multiplying through by a factor which makes £,,; = 1. 
If this last feature is incorporated into a program, the whole computation can be 
expressed as the following iterative algorithm: 


A= batt, 161 + bn41,2&2 +++ + Dn+1, nén + Ontt, n¢1 


1 
ti! = d {brits + Dioge +--+ + Ding + 51, n41} 
1 
(1) fx! = 5 {baits + bake + -+- + danke + ba nts} 
1 
é, 7. r {barks + Dnoks o> + Dankn + Bn, nea}. 


The algorithm is repeated until each é,’ agrees with &; to the required degree of 
accuracy. The last \ and col (£1, 2, --- &,, 1) are the required proper value and 
vector respectively. 

This algorithm has appeared in numerous places, e.g., von Mises and Pollaczek- 
Geiringer [1]. 

The program thus set up can be directly applied to a system of a linear equa- 
tions in unknowns as follows. The system of equations AX = b, where A isa 
non-singular (m,) matrix with real entries and X and 6 are the real column 


vectors col (X1, Xo, ---, Xn), col (b1, be, ---, bn), respectively, is equivalent to 
the system 
A+kI —b al 
0 , k y 7 'y 
where y is the vector col (X1,; X2, ---, Xn, 1) and & is any real constant. This 


shows that y is a proper vector of the matrix 


_(A+kl -b 
B-( 0° =" 


corresponding to the proper value k. 

The iteration corresponding to the equations (1) will converge to y pro- 
vided || is larger than the modulus of any other proper value of B. The existence 
or non-existence of an appropriate k is covered by the following theorem. 

THEOREM 1. Let A be non-singular. A necessary and sufficient condition that 
there exist a real number k such that each eigenvalue distinct from k of the matrix 


_(A+kI -b 
B-( 0° ite 
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be of modulus less than |k| is that all eigenvalues of the matrix A have real parts 
different from zero of the same sign. 

Proof: Necessity: Let A have proper values Ax, A2, ---, An. Then B has proper 
values k, k + Ax, & + Ao, ---, & + An. If two of the Aj, e.g., Ax, A2 have real parts 
of opposite sign, then |& + A:| <|&| implies that the real part of \ has sign oppo- 
site to that of k. Hence )2 has a real part of the same sign as that of & so that 
|k + A2| >|]. Also, if any A; has real part 0, then |k + A;| >|%| unless |A| = 0, 
for all real k. 

Sufficiency: If X1, 2, ---, An have real parts all of the same sign, first suppose 
this sign to be negative. In the complex plane (see figure 1) each X, is to the left 
of the y-axis. Let u; be the point at which the line joining the origin to A; meets 
the circle |z + 1|= 1. Let k be any real number such that k >|A,|/|u;| for 
4 = 1,2, ---,m. Then ,/k lies inside the circle |z + 1| = 1. Hence |1 + A,/k| < 1, 
ie., |& + ;| < & for all 7. In the case where the real parts of \; are all positive 
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pa 








Fic. 1 


it is only necessary to replace the circle |z + 1| = 1 by the circle |z — 1] = 1 
and then choose k in such a way that — k >|A;|/|u:| for all 2. 

In general, while criteria exist for deciding whether the real parts of the proper 
values of a matrix A are all of the same sign, such criteria are not computationally 
of much use. If it is known that the proper values of A are all real (e.g., in the 
case where A is symmetric), then if a value of & exists for which the iteration 
converges, the value k = — tr A will be adequate. In any case the sign of a 
suitable & will always be opposite to that of tr A. 

From a practical point of view for any A a simple procedure is to assume a 
value of k exists and choose a large multiple of — tr A to be this value. If the 
iteration does not appear to converge, it is suggested that the equations be 
modified by the usual type of linear transformation. If one cannot find an obvious 
row transformation which makes the system converge it may be noted that, if A 
is non-singular, A7A is positive definite so that the system AX = 6 may be 
replaced by A7AX = A‘). For this system k = — max a,;, where a;; are the 
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elements of A7A, is always adequate. While as Taussky shows [2], the system 
ATAX = A’ is more ill-conditioned than AX = b, this ill condition reflects 
itself in a larger computing time but does not affect accuracy. 

The system has been applied with success at the Canadian Armament Re- 
search and Development Establishment on the ALWAC III machine. The cir- 
cumstance was such that the routine for proper vectors was programmed prior to 
that for matrix inversion. It was originally meant as a stop-gap measure but has 
turned out to be a practical method of much merit. 

N. S. MENDELSOHN 
Canadian Armament Res. and Dev. Est. 
and 
The University of Manitoba 
Winnipeg, Canada 


1. R. von Mises & H. PoLtaczeEK-GerrRincer, “Practische Verfahren der Gleichungs- 
auflésung,”’ Zeit. angew. Math. Mech., v. 9, 1929, p. 58-79 and 152-164. 

2. OtGa Taussky, “‘Notes on numerical analysis—2, Note on the condition of matrices,”’ 
MTAC, v. 4, 1950, p. 111-112. 


The Computation of Complex Proper Values and Vectors of a 
Real Matrix with Application to Polynomials 


Introduction. The power method for computing the largest real proper value 
and corresponding proper vector of a real matrix A can with very minor adjust- 
ments be applied to the complex proper values and vectors. A good exposition of 
the power method may be found in von Mises and Pollaczek-Geiringer [4]. 
When applied to the companion matrix of a real polynomial the method yields 
a real root or pair of complex roots of the polynomial. It turns out that in this 
application the method is equivalent to Bernoulli’s method, as was shown by 
Aitken [1]. Having obtained the root or roots of maximum modulus other roots 
are usually obtainable by the elementary devices indicated below. The method 
admits of easy programming on a digital computer. The program has been carried 
out on the ALWAC III computer at the Canadian Armament Research and 
Development Establishment and has turned out to be very practical indeed. 

Method and Theory. Let A be an (r, 7) real matrix with proper values A, 
he, «++, Ay and corresponding column proper vectors x1, x2, ---, x,. It is assumed 
that |Ai| > |A2| > |As| > --- > |A-|. If any one of the A; is not real its conjugate 
appears among the remaining A; and the corresponding proper vectors are 
conjugate. If » is an arbitrary non-zero column vector with r components, 
n can be expressed in the form 9 = cx; + Cov. + --- + ¢,x,. It follows that 
Qn = A™ = C1A"X1 + Cod2"X2 + yg + Crds"Xy. 

If |Ai| > |As|, then A; is real and 


= ax fe +a(™)"s +o +e(*) af. 
Nn = Al 1%1 2 i 2 r i r 


If c, + 0, the dominant term in 7, is A,"¢:x,. Hence, for sufficiently large n, 
‘nti = Asma. It is ordinarily convenient at each stage to normalize 7, so that the 
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last component is 1. Accordingly putting 7, = col (y{", yf” --- y/®,, 1), the re- 
lationship between 7, and 7,4: is given by: 


) 
Ba+i = any” + ays” + S49 + Dr, 1 ont + Orr 


1 
yt) Te fasys” + ays” +03 + 01, 7-191 + Gir} 
Bn+ti 








1 n n n 
yt? = {aaxry} ' + axy$ ; +-°° + G2, r1y + a2} 
(I) Mn+1 
yt? = {dy—1, wy” + ar-1, ws” + --- + @y-1, 1s + @,—1,r}. 


Mn+i1 


Here a,; is the (i, 7) element of the matrix A, un»; is a factor of proportionality 
which makes the last component of 9,4: = 1, and the iteration is started using a 
non-exceptional real vector 7 = 0 with last component 1. The vector 7, converges 
to the proper vector x; and un»+1 converges to A1. 

In the case where |Ai| = |A2| >|As|, A1 # Ae the iteration (J) no longer con- 
verges. This case occurs when either A; = — A2 or when A; and A: form a conjugate 
complex pair. We will assume the latter occurs. Put A2 = \, and x2 = 2,. Although 
the iteration (J) diverges it is still useful, as is shown by the following theorem. 

THEOREM. Let A be an (r,r) matrix with real entries and proper values 1, 
Xi, As, Aa, «+ * Ay Gnd corresponding proper vectors x1, £1, x3, ---,X,- such that d, is 
not real and |i|=|Ai| > As D> Ag > --- D Ay. Let unys be given by the iteration (I), 
where (I) ts started by means of a real, non-zero vector no which has 1 as its last 
component. Then the sequence 


(Un — Mn+ 2) 


Vn = $unt1 
(Un - Mn+1) 


converges to a limit a. Furthermore, the sequence 


(unt. — Mn+2) 


on = Mnbn+1 
(tn — Mn+1) 


converges to a limit p* and the proper values dx, \1 are given by a + i8 where 
B = Vp? — a’. ‘ 
Proof. Let no = ¢1%1 + Cov, + --- + .¢,x,. Then 


tn = Rn{cidr™%1 + Code™X%2 + +++ + CAr"x} 


where , is a constant which makes the last component of 7, = 1. For sufficiently 
large m all but the first two terms of the sum may be neglected. Let N be a specific 
value of m sufficiently large so that yy = ax; + G%, the neglected terms being 
negligible. Since the vector x; is only determined up to a constant factor we may 
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replace x; by ax;. With this convention ny = x; + 2:. Let £, be the last component 
of x:. Since ny has last component 1, £ + £ = 1. By the iteration (J) it follows 
that nw41 = Ann/un+1 = (Axi + Ai1)/uw41. Also, since 94: has last component 
1, uaa = Aake + AE, = (Anke + Ank,)/(E- + E). Put Ar = a + 18 = pe*, t, = der, 
cos (26 + y) 
cos (6 + 7) 
On eliminating y the relation p? — 2uy+1p cos @ + uwiiuwi2 = O is obtained. Re- 
placing N by N +1 yields p* — 2un429 cos 0 + uns2unis = 0. Solving these 
latter two equations, one obtains 


Then py+1 = pcos (6 + +). In the same way it follows that uvs2 = p 


. Lpewse(unsi — wv4s) 
a = pcos@ = — = yyi1; p’ 
2 MN+1 — BN+2 MN+1 — BN+2 





= My+isn+2(un+e oni Hn+3) 





ON+1- 


To solve the problem of determining A; on a digital computer, the numbers 
zn are first generated from the iteration (J) and then the », and o, are computed 
from the yu». When successive values of v, and o, differ by less than the allowable 
preassigned error the iterations are stopped and a and @ are given the values 
a = Yny1, B = Vony1 — Yex1. This algorithm for obtaining a, 8, and p is also dis- 
cussed in Rutishauser [3] and Aitken [1 ]. 

A word of warning may be appropriate here. While theoretically the », and 
gn converge, it may happen in practice that rounding off errors in computed 
values of these will leave the magnitude of the difference of successive values 
greater than the preassigned error. Hence, if the computation continues for too 
long a time it may be necessary to print out a few consecutive values of », and 
o, and examine the differences. 

In the case of complex proper values, there may be no occasion to ob- 
tain the corresponding proper vectors, but the computation presents no diffi- 
culty. Having computed A, the vector x, may be obtained from the equation 
x, = k(unsimw+1 — Arma) where & is a constant of proportionality. 

Application to Polynomials. Let f(x) = x* + cx" + cox”? + --- +6,=0 








have roots Ax, Ao, «++, Ae such that |A,| > /A.|---2>]A,-|. The companion matrix 
f 0 1 0 0 0 ) 
0 0 1 0 0 
Aa 
0 0 0----0 i 
* —Cry —Cr—1, —Cr—2, > Gh 
has proper values Ai, A2, -:-, Ar. The computation method here described will 


yield either a single real root or a pair of complex roots. Further roots of f(x) = 0 
may be obtained by various devices. For instance the matrix A + kJ has proper 
values A; + Rk, Ao + Rk, «++ Ae + &. It may well happen that |A; + &| is not the 
largest of the numbers |A: + &|, |A2 + 2], ---, |Ay + &]. In this case, another 
root or pair of roots of f(x) = 0 is obtained. If A; is real the choice k = — A, is 
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certain to yield further roots. If 4; = \, are a pair of conjugate complex roots, 
taking k = — real part of \; will minimize |\, + | and will possibly yield further 


1 
roots. By forming the polynomial g(x) = x’f ( ) one obtains the roots of mini- 


mum modulus of f(x). These devices have the advantage that the precision ob- 
tained in any new root is independent of the errors in previously obtained roots. 
It is only after these devices have become exhausted, that it may become neces- 
sary to remove from f(x) factors corresponding to computed roots. 

Concluding Remarks. The method described here fails only in the case where 
f(x) has more than one pair of roots of the same maximum modulus (e.g., the 
equation x” — 1 = 0). It could be refined without difficulty to take in these more 
general cases, as indeed was done by Rutishauser [3] and Aitken [1 ]. The results 
are considerably more complicated in form and it is doubtful if they have any 
practical value. 

N. S. MENDELSOHN 
Canadian Armament Res. and Dev. Est. 
and 
The University of Manitoba 
“a Canada 
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H. RUTISHAUSER, “Anwendungen des Quotienten-Differenzen- Algorithmus,” Zeit. angew. 
Maih. Physik, v. 5, 1954, p. 496-508. 


4. R. von MIsEs & H. PoLvLaczeEK-GEIRINGER, ‘“‘Practische Verfahren der Gleichungs- 
auflésung,” Zeit. angew. Math. Mech., v. 9, 1929, p. 58-79 and 152-164. 


A Method of Inverting Large Matrices of Special Form 


1. Introduction. A method is suggested to invert large matrices, M,, of the 
form 


ao —P, 0 0 vee 0 
—P, a —P, 0 "es 0 
ae 0 —P, ae —P; 
(1) M, rw 0 0 —P; a3 
—P,, 
0 0 —P, an 








in which the elements are themselves special square matrices. The diagonal 
elements, a;, are symmetric matrices. The minus signs are not required; they 
were a convenience in the particular problems we studied. This type of matrix 
arises, e.g., in the least-square fitting of survey data or in the difference equation 
approximation, to some common partial differential equations. 

The method suggested takes advantage of the special properties of the large 
matrix and involves only algebraic operations and inversions of matrices the size 
of the elements. It is thus particularly applicable for use by ‘‘medium”’ sized 
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computers. The routine will be derived using the simpler matrix, 


95 





a -I 0 0 
—I ai —I 
0 —I ae —I 
(2) M, = 
0 
—I 
0 0 —-I a 








where J is the identity matrix. 

This arises in the difference equation solution to Poisson’s equation over a 
rectangle as shown by Stein and Peck [3]. The algorithm is then generalized_to 
include the case where P; are diagonal matrices. It could be further generalized 
to include P; with no special properties if desired. 

A method of partitioning is used. The basic algebra of the method is presented 
in Frazer, Duncan, and Collar [2]. The method of partitioning by one row and 
column at a time was suggested by Wagner [4]. A comparison of partitioning 
with other methods of matrix inversion is given by Fox [1], who also gives a 
comprehensive bibliography. 


2. Method of Inversion. Let 


: u- (at) 


where a and d are square matrices. Let the inverse of M be given by the equation 


AIC 
_ U . 
(4) - (415) 








Then according to H. M. Wagner [4], 


D = (d — ca)" 
C = —a"bD 
B= — Dea 
A =a" —a"bdB. 


Now define the set of sub-matrices 
= ao —I me a, b; 
©) WE oe a ) ee (° rd 


( " 

—I al > *) 
= 0 —I ae - Ce d, , 
and, in general, 


( Mi. | 0 ) 
—I a; rt 
(6) ms = 0 =f ai 2 (* d; ; 
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Application of equations (5) obtains 


Mt ( + ag? (ay — ao") ag ag (a, — ws 
r= , 
(a: — ao™)tag (a; — ag) 


This suggests a convenient notation 


[0] =a? 
[1] = @ — (0) 
(7) ' 


[n] = (a, — [m — 1])-. 
In this notation, 


_, _ ((01+ (oIE13C0] c07C17) 
7 ( [1 [0] 4 


Using this expression and applying equations (5) to M2 gives 


CO} + COJC1 ILO} + COILIIL2ILAILO}, (COWL1] + COWII2I01), COILI2) 
(8) Mt = ( Cj00]+  Mj2I1j(0], C1J+ (Ctj(2I01], (11023). 
(2I(1(0], (27C1], (2) 


Using mathematical induction and equation (8) as a guide, it is easy to write 
down M,,- in terms of the matrices [7] for any value of m. An explicit expression 
for a general term in the inverse is not apparent. The following recurrence rela- 
tions proved convenient for machine programming. Let the elements of M, be 
m;;(4,j7 = 0,1,2--- m). 


Then 
Mnn = [n] 
(9) mis = [t]mini;, t> 7 
ni = [t] + My, i441 [i] 
m,; = transpose m,;, 14 < j. 


3. Generalization of Result. By similar analysis, an expression for the inverse 
of the matrix of equation (1) may be derived. The derivation itself is not of 
interest. For the general case in which the matrices P; are general, equations 
(7) are replaced by 


[0] = ac 


[1] = (a: — P,L0JP:)* 
(7’) 


[n] = (an — Paln — 1]P,)-. 
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Equation (8) becomes (omitting brackets) 


0+0P,;1P,0+0P)1P.2P.1P,0, 0P\1+0P,\1P22P.1, 0P,i1P.2 
(8’) M;" = 1P,0+ 1P.2P.1P,0, 1+ 1P,2P21, 1P,2 
2P21P,0, 2P,1, 2 


Recurrence formulas similar to (9) may be induced from this array. 


4. Computation Procedure. We have inverted matrices of both described types 
using a Datatron computer without magnetic tape storage. It was found con- 
venient to use interpretive matrix algebra commands of a three-address form. 
Storage locations were reserved both for square and for diagonal matrices. These 
locations were given pseudo-addresses. Typical single word commands would 
perform the following operations: 


(i) Multiply the matrix in pseudo-address A by that in B and store the 
product in C 
(ii) Invert the matrix in A and store the inverse in B 
(iii) Punch onto cards the matrix (or its transpose) stored in A. 


A standard punch card format is used. 

The sub-matrices [0] to [#] are first computed and punched. These results 
are used by a second program to compute the inverse. Using the interpretive 
three-address commands, a large inversion may be programmed in a few minutes. 
Our programs limit us to:20 X 20 matrices (a; and P;) and to a 200 X 200 ma- 
trix M. It required about 2} hours to invert an 80 X 80 matrix by this method 
on our computer. 

H. I. MEYER 


B. J. HOLLINGSWORTH 
United Gas Corporation 


Shreveport, Louisiana 
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3. P. Stern & E. L. PECK, “On the numerical solution of Poisson’s equation over a rectangle,” 
Pac. J. of Math., v. 5, 1955, p. 999-1011. 

4. H. M. WAGNER, Matrix Inversion on an Automatic Calculator by Row and Column Partition- 
ing, RAND Corp. Report P-417, 1953. 


Note on the General Solution of the Confluent Hypergeometric Equation 


The note [1] on a Webb & Airey-Adams-Bateman-Olsson error by Murlan 
S. Corrington, which is concerned with the confluent hypergeometric differential 
equation 


(1) xy” + (y — x)y’ — ay = 0, 


does not quite cover all the possibilities for its general solution. A statement of 
the complete situation does not seem to have appeared in print before now, though 
in the Fletcher, Miller, and Rosenhead Index of Mathematical Tables [2], the 
statement is deliberately cautious and correct. It seems worth while to give a 
full statement here. 
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If a and ¥ are constants, which may be complex, the complete solution of (1) is 
y = AM(a, 7, x) + Bx! 7*M(a — y + 1,2 — 7, x) 


where A and B are arbitrary constants, and 


Néine wripts g fOtDZ , ... 


yi! yv(y +1) 2! 
so long as both solutions exist and differ. 

The first solution ceases to have a meaning when y¥ is zero or a negative in- 
teger, unless a is also a non-positive integer such that |a| < |y|. In the latter case 
the numerator vanishes with or before the denominator, and both solutions still 
hold. In other cases, as y approaches zero or a negative integer, M(a, y, x)/T (y) 
remains finite, but becomes a multiple of x'-7M(a -—- y + 1,2 — y, x). 





The case a = y = — p, p zero or a positive integer, as indicated in the Index, 
is of peculiar interest. The complete solution then becomes 
A. 008 xP xt apt? ) 
yoA(ttere 4% 4 ‘ += )+3( 4+ @+a!i* da 
x? 
= (A ~B(ite+%4 vee +=) + Be. 


Since the substitution y = x'-7z interchanges the roles of the “‘first’’ and 
“‘second”’ solutions, leaving a differential equation that is still confluent hyper- 
geometric, and since one of y and 2 — vy, here assumed real, must be positive, 
it is clear that there is no loss of generality if y is taken positive; the logarithmic 
form of solution thus need not be considered for the‘first solution. 

The second solution is not distinct from the first if y = 1 and ceases to have 
a meaning if y 2 1 is an integer, unless a is also an integer such that 1 < a < ¥. 
In the latter case, as before, the numerator vanishes with or before the denomi- 
nator and both solutions are valid. 

The cases emphasized by use of italics, often overlooked (as by Corrington) 
lead to polynomial M-functions. 

When y 2 1 is an integer, whilst a is not a positive integer less than y, the 
second solution of (1) is replaced by a logarithmic solution, given by equation 
(2) of Corrington’s note; but this is not the complete solution, as stated, a further 
term AM (a, y, x) being needed. This addition, with A replaced by a new arbitrary 
constant D = A + C{y(i — a) — ¥(y) — ¥(1)}, and with the Beta-function 
terms written out, replaces Corrington’s solution by the full solution 


y = (CInx + D)M(a, y, x) + CN(a, vy, x) + CS(a, y, x) 
in which 
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and 
S(a, 1,2) = (-1)T@ — 7 + Py — 1) 2 a 
I' (a) 
a-yti«x (a-yt+i)@—-—v+2)# ~— 
x (14 2-7 1! (2 — y)(3 — y) at 


toy — 1 terms) 
y- 11) G-DO-2)1, @-1)GO-9)G-392!_ 
a-1ix (a—i1)(@—2)x (@—1)(a — 2)(a — 3) x 


(y—1)(y—2)++-24  (y— 2)! 
(@@—1)(@-2)---@-y+1) 








+ (—1)7 


This function S(@,7,x) is the part omitted by Webb and Airey, and others. 
Of course, if y = 1, then S = 0 and the solution given by Webb and Airey is 
correct. 

J. C. P. MILLER 
The University Mathematical Laboratory 
Cambridge, England 


1. A Webb & Airey-Adams-Bateman-Olsson Error, M.T.E. 117, MTAC, v. 2, 1947, p. 352-353. 
2. A. FLetcuer, J. C. P. Miter, & L. RosENHEAD, An Index of Mathematical Tables, Scien- 
tific Computing Service Ltd., London, 1946, article 22.56, p. 336. 


Approximation and Table of the Weierstrass g Function in 
the Equianharmonic Case for Real Argument 


The most readily accessible table of the Weierstrass g function for the equi- 
anharmonic case appears to be that in Jahnke-Emde [2], which is mainly to 4D. 

The equianharmonic case for the function g(u; ge, gs) is that in which g2 = 0, 
gs = 1. The function behaves like 1/u* near the origin. This leads us to consider 
tabulating 


(1) f(u; 0, 1) = g(u; 0,1) — 1/2. 
Because of symmetries, we may restrict attention to the interval 0 < u < ws, 
where w: ~ 1.52995 and is the real half-period of g; see F. Tricomi [1]. 


It is possible to represent g on the above interval to 7S by a relatively simple 
approximant. Write 


(2) f(u; 0,1) = cgut + cou! + cou®® + Creu + crgu*®, 


the polynomial being obtained by truncation of Maclaurin’s series. Thus we have 


(3) e(u; 0,1) = g(u)/u* = h(y)/*, 
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where g(u) = 1 + wf(u; 0,1), y = u®, and h(y) is of the fifth degree. It turns out 
that hk can be ‘‘economized by the use of Chebyshev polynomials” in the sense 
of C. Lanczos [3] so that 7S values of g can be obtained on the interval of 
interest using a cubic approximation to h(y). That is 


(4) @(u;0,1) ~ H(y)/w, where y = uw’ and H(y) = = .99999 990 
+ .03571 4504y 
+ .00009 80386 9? 
+ .00000 01932 Oy’. 


At the half-period u = w2 ~ 1.52995, H(y)/u? ~ 0.62996 055 whereas the true 
value for g(w2) = 4-* ~ 0.62996 052. 

A compact printed table of the auxiliary function f has been constructed for 
use with Everett’s Interpolation Formula to second (modified) differences. It 
follows. The differences have been modified by throwback of fourth differences 


Table of Weierstrass Function f = g — 1/u’. 


u f(u; 9, 1) 6*f u f(u; 9, 1) o*f 
0.0 .00000 00 - 85 0.80 .01463 91 68 89 
0.1 .00000 36 + 341 0.85 .01866 24 77 95 
0.2 .00005 71 16 30 0.90 .02346 64 87 66 
0.3 .00028 93 37 70 0.95 .02914 83 98 10 
0.4 .00091 43 67 71 1.00 .03581 26 109 27 
0.5 .00223 22 106 32 1.05 .04357 12 121 30 
0.6 .00462 92 153 53 1.10 .05254 46 134 32 
0.7 .00857 78 209 57 1.15 .06286 32 148 39 
0.8 .01463 91 274 77 1.20 . .07466 81 163 77 
1.25 .08811 35 180 64 
1.30 -10336 86 199 30 
1.35 .12062 07 220 11 
1.40 .14007 87 243 53 
1.45 .16197 78 270 08 
1.50 .18658 49 300 52 
1.55 .21420 60 335 70 


using the throwback constant —0.18393. This throwback process is described by 
various authors; see for example, the NBS work [4] or F. B. Hildebrand [5]. 
The formula for H(y) has been published previously without comment (see 
Review 58, p. 110 this issue MTAC). 
Tuomas H. SouTHARD 


University of California 
Los Angeles, California 


This paper was prepared under the sponsorship of the Office of Naval Research. 

1. F. Tricomi, Elliptische Funktionen, Akademische Verlagsgesellschaft, Leipzig, 1948, p. 20. 

2. E. JAHNKE & F. Empe, Tables of Functions with Formulae and Curves, 4th ed., Dover 
Publications, New York, 1945, p. 100. 

3. NBS Applied Mathematics Series No. 9, Tables of Chebyshev Polynomials S,(x) and C(x), 
U. S. Gov. Printing Office, Washington, D. C., 1952, p. 15. 

4. NBSCL, Tables of Bessel Functions of Fractional Order, v. 1, Columbia University Press, 
New York, 1948, p. 33. 

5. F. B. HILDEBRAND, Introduction to Numerical Analysis, McGraw-Hill Book Co., Inc., 
New York, 1956, p. 112. 
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A Lemma for Automatic Optimum Programming on the IBM 650 


1. Introduction. The important problem of optimum programming a set of 
instructions for the IBM 650 Magnetic Drum Data Processing Machine has been 
attacked in several ways [1, 2, 3, 4]. The techniques commonly used, whether 
they be manual or machine operated, refer to the optimum programming table 
for the minimum spacing of instructions and data on the drum [5]. Perhaps the 
major limitation in current procedures is that the program becomes less nearly 
optimum as succeeding instructions are placed; i.e., unless the program contains 
only relatively few instructions that need not be densely located on the drum, 
once half of the program has been placed optimally the second half is located on 
a “catch-as-catch-can” basis and reference to the optimum programming table 
for placement becomes virtually useless. ‘ 

It is the purpose of this paper to present a lemma which has proved useful 
for automatic optimization of 650 routines. 

2. Fundamentals of 650 Coding. The main storage medium of the 650 is a 
revolving (12,500 r.p.m.) magnetic drum containing 1000 or 2000 locations, which 
hold signed 10-digit words of data and instructions. Drum locations, assigned 
addresses 0000, 0001, ---, 1999 (assuming a 2000 word drum), are arranged in 
sequential order, starting at 0000, in strips of 50 words around the surface of the 
drum. In addition to drum storage, there are four other special locations, address- 
able as 8000, 8001, 8002, and 8003. 

A 650 instruction is composed of three parts, the two most significant digits 
specifying the operation, the next four digits indicating the Data Address, and 
the last four digits indicating the next Instruction Address (the sign of an in- 
struction has no effect on the operation). The Data Address usually denotes the 
location for a) obtaining an operand, b) storing a result, or c) reading or punching 
information; the Data Address may also specify the length of a shift or the address 
of the next instruction that follows a logical test in which the test condition (e.g., 
non-zero accumulator) has been satisfied. The Instruction Address usually refers 
to the location of the next instruction, the exception being after a logical test. 
Both the Data and Instruction Addresses may be any drum location or 800X, 
with the qualification that 800X may not be used as a Data Address in input- 
output, store, or table look-up commands. 

The presence of a next Instruction Address allows one to disperse instructions 
as well as data around the drum; hence for the most part, coding a sequence of 
required operations may be undertaken without regard to the ultimate placement 
of instructions and data. Storage is finally assigned to reduce waiting time as 
much as possible, i.e., the time between the moment the computer is ready to 
receive some data or an instruction and the moment the data or instruction arrives 
under the drum read-write heads. Because the sequence of operations and the 
assignment of storage are practically distinct problems, minimum latency coding 
in the 650 is conceptually easier to handle than, say, in mercury delay line sys- 
tems, where instructions are placed in successive locations, and, as a result, the 
sequence and assignment questions become inseparable. 
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3. The Optimizing Lemma. Let and k be positive integers. Then the sequence 


a) 1,1 +,1-+ 2n,1 + 3n, --- modulo (nk + 1) completely exhausts the 
numbers (1, 2,3, ---,#k +1) before repeating any number in the sequence; 
similarly for 

b) 1, 1+, 1+ 2n, 1+ 3n, --- modulo (m(2k + 1) + 2)/2 where a is 
even; and 

c) 1,1 +, 1+ 2n, 1 + 3n, --- modulo (m(2k + 1) + 1)/2 where 1 is odd. 


For example, let nk + 1 = 10 and m = k = 3; then we have the sequence 
1, 4, 7, 10, 3, 6, 9, 2, 5, 8, which completely exhausts the numbers 1 through 10 
without any repetition. 

If, in an effort to reduce access time on the 650, we desire to place a sequence 
of N consecutive instructions and data (a precise definition of “consecutive’’ is 
given below) m locations apart, we find an integer k such that one of the cases in 
our lemma applies. Thus we may place a set of 50 consecutive instructions and 
data 7 locations apart on a band of 50 locations: 0000, 0007, 0014, ---, 0049, 
0006, ---, 0043 (n = k = 7, case a) of the lemma). Furthermore we may place 
on the drum any set of N consecutive instructions and data by first densely filling 
up r* bands of 50 locations, where r* is the largest integral value of r such that 
N — 50r > 0. The remaining N — r*50 instructions are placed according to the 
lemma within the next block of 7k* + 1 successive locations, where k* is the 
smallest integral value of k such that 7k + 1 > N — 50r*. If storage space is not 
limited a better choice is k* = 7. As the reader can see, the previous example 
generalizes to any N and and need not utilize integral multiples of 50 locations. 

4. Sequential (or Serial) Coding. We shall now clarify the notion of a 
“sequence of consecutive instructions and data.”’ Consider a given 650 instruction 
stored in location j. If the instruction’s operation*is add, subtract, multiply, 
divide, store, or load, we shall code the Data Address as location 7 + 1 and the 
next Instruction Address as location j + 2. If the Data Address is 800X, we may 
code the Instruction Address as j + 1, or as 7 + 2 and void the location j + 1; 
of course the Instruction Address may be 800X if desired. On branching operations 
we should let location 7 + 1 be the most traversed route; the address remaining 
will refer to some other part of the program. In shift operations we code the 
Instruction Address as location 7 + 1. Most routines will not have all their in- 
structions and data in strictly consecutive order, since branching operations may 
be included and since some data locations may be referred to more than once in 
a program. The optimizing routine suggested here operates best if the program 
is originally written as nearly sequentially as possible. 

The coder observes the following procedure in writing a new routine: He codes 
his program in as nearly a “sequential order” as is practicable with real drum 
addresses starting, say, at s; and ending at s,, and has it punched on any standard 
one/card form. He next debugs the sequential code on the 650. As experienced 650 
programmers will attest, the sequential pattern greatly aids the checking out 
procedure. Any additional operations which must be included are added in loca- 
tions immediately preceding s; or following s.. The final debugged program is then 
located between, say, sy and s,; note that the program need not be dense between 
sy and s,. The one/card form contains in addition to the location of the instruction 
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and the instruction itself, coded information analogous to that needed for auto- 
matic translating routines, viz.,a number XX (99, 98, 89, or 88) referring to 
whether or not the Data or Instruction Addresses are translatable (9 being 
not-translatable). 

5. Automatic Optimizing. To program the debugged routine, the coder inserts 
in the optimal programming deck an “information card’’ which has punched in 
it Sy, Sg, and 0, (the first location for the optimal program). He places the one/card 
debugged routine behind the optimizing deck and runs the combined deck through 
the 650; a new one/card optimized deck is punched out. 

The optimizing routine actually consists of two separate parts. 

We shall assume that some m and some modulus have been selected to satisfy 
a particular case in our lemma. In the first part, given the addresses s,; and sy, 
the 650 computes how many locations are needed for the optimized routine. 
Starting with address 0,;, the machine constructs the sequence of addresses as 
prescribed by our lemma, and concomitantly, the addresses in the new sequence 
are associated with and correspondingly stored in locations s, through s,. When 
the first part of the program has been completed, the second part commences by 
reading each card in the one/card debugged deck, optimizing its contents, and 
subsequently punching the new optimized one/card instruction. From the 9’s 
and 8’s information punched in the original card, each instruction is altered so as 
to replace where indicated a Data or Instruction Address in the range sy; to s, 
with its associated new optimal address. The correct new location of each instruc- 
tion in the sequential deck is determined by testing whether the old instruction 
location lies within the range sy, to s,; if it does not, the location is not altered. 
Hence cards used in the sequential program deck with instructions and data to 
be placed in absolute (i.e., fixed or unalterable) addresses may remain in the deck 
which is being optimized ; although the fixed address will not be altered, necessary 
Data and Instruction Address changes are made as indicated by the 9’s and 
8’s coding. 

We may now summarize the programming limitations imposed by this 
method of optimizing: No absolute address may lie either within the range sy, to 
S, OF 0s to o, (the last optimized location). Depending on the particular optimizing 
program devised, the range sy; to s, may or may not be allowed to overlap with 
the locations in which the optimizing program is stored. There are no other 
restrictions on Sy, Sy, 0s, 0,, other than that they must be legitimate 650 addresses, 
and s; to s, may overlap 0; to 0,. 

6. General Comments. It should be obvious to the reader that the proposed 
scheme has the property that the degree of optimization is essentially invariant 
to the number of instructions previously optimally located. Thus in a program 
utilizing 500 locations, the instructions and data can be placed optimally and 
densely in 500 consecutive locations, and the degree of optimization for the last 
100 instructions and data is essentially the same as that for the previous 400 
instructions and data. It is the author’s belief that this property will prove to 
yield results at least as good as those obtained by ‘‘strict’’ optimization of the 
instructions in the first half of a program and ‘‘catch-as-catch-can” allocation of 
the instructions in the second half. 

Using D. W. Sweeney’s statistics [6] for the average use of 650 commands, 
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we may expect our lemma for m = k = 7 (case a)) to be about three times as 
fast as sequential coding. (In one program coded by the author, the optimizing 
of a payroll calculation resulted in a reduction in computing time sufficient to 
cause calculations to take place at full reading and punching speeds.) A desirable 
property of the suggested optimizing code is that a single routine may be written 
for one case of the lemma and arbitrary m and modulus; assignment of particular 
values to these parameters may be done by the insertion of a single card in the 
optimizing deck. Therefore a sequential program may easily be optimized for 
several different values of m and modulus in order to test the best mode of optimiz- 
ing for the particular routine. If a program contains mostly add, subtract, store, 
and branch commands, computing speed may be increased by using m < 7. 
Given 1, one should select a case of the lemma and k to maximize (modulus) < 50. 

Although the lemma is simple enough to use as a manual technique for con- 
comitantly coding and optimizing, we do not recommend the procedure, because 
testing out a sequential routine has important time saving advantages and because 
the automatic method for optimizing is so easy and fast. 

A routine using case a) of the lemma has been written by the author for 
M.I.T.’s Statistical Services; the program occupies less than 100 locations and 
will optimize a one/card deck at full punching speed of 100 cards/minute. 


HarvEY M. WAGNER 
Massachusetts Institute of Technology 
Cambridge, Massachusetts 


1. B. Gorpon & A. J. Datton, “Optimizing program,” Equitable Life Assurance Society 
and IBM, New York, 1955 (mimeographed). 

2. Barry Gorpon, “An optimizing program for the IBM 650,” J. Assn. for Comp. Machinery, 
v. 3, 1956, p. 3-5. 

= = POLEY & G. L. Mitcuett, “S. O. A. P., IBM 650 symbolic optimal assembly program, 
programmer’s guide,” IBM, New York, 1955. 

4. E. F. SHEPHERD, “An automatic method of optimum programming,’’ JBM Technical 
Newsletter No. 10, New York, 1955, p. 95-104. 

5. IBM, “Type 650 manual of operation, first revision,’’ New York, 1955, p. 70-89. 

6. D. W. SWEENEY, “A note on optimum programming and the IBM type 650 operation code 
usage,” IBM Technical Newsletter No. 10, New York, 1955, p. 105-107. 
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49.—A. V. LeBEDEV & R. M. FEporova, Spravochnik po Matematicheskim 

Tabliisam (Guide to Mathematical Tables), Moscow, 1956, xlvi + 549 p., 

26 cm. Price 29.20 rubles. 

In this Guide to Mathematical Tables the authors have extended the largest 
available index to mathematical tables, Fletcher, Miller, and Rosenhead [1], 
to include tables published in books up through 1952 or later and tables published 
in journals through 1953. These dates seem more recent than the closing dates 
for Schiitte’s index [2]. 

The preparation of this index was necessitated by the great activity in table 
making and related parts of numerical analysis in USSR. The Guide was prepared 
as a convenience both to scholars and to practical workers. Thus, it is to serve a 
somewhat wider class of users than Fletcher, Miller, and Rosenhead [1] or 
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Schiitte [2]. On the whole, [1] was directed to the fairly advanced computer, 
and references to elementary tables other than the most basic were largely 
omitted. For users outside USSR, the present Guide will serve approximately the 
same purpose, but four- and five-digit tables most easily available in USSR were 
listed and some items not conveniently available in USSR (or elsewhere) are 
omitted. The present work seems to share with Fletcher, Miller, and Rosenhead 
[1] the decision to eliminate tables based on physical measurements. Schiitte 
[2] directed his effort toward the listing of practically useful and available 
tables. 

There are presumably no entries in Fletcher, Miller, and Rosenhead [1 ] which 
are accidentally excluded from the present Guide, for the authors acknowledge 
having used [1], and also point out that many of their references are from other 
indexes or references and are not from a direct examination of the table reported. 
They also acknowledge help received by consulting H. T. Davis [3] and MTAC. 
However, many references in [1] are not found here—five of the twenty-five 
references to J. R. Airy noted in [1] are not found here. Presumably this is in 
accord with restrictions imposed purposely by the present authors. 

A detailed table of contents extending over 42 pages makes reference to the 
book somewhat easier than reference to earlier indexes. However, the casual user, 
who has not memorized the classification system, will lose much time leafing 
through this table of contents seeking proper chapter headings. There is a name 
index, but no index of functions tabulated. 

Approximate translation of the chapter headings follows: 


Chapter 1. Powers, and rational and algebraic functions 

Chapter 2. Trigonometric functions and various quantities connected with 
circles and spheres 

Chapter 3. Exponential and hyperbolic functions 

Chapter 4. Common and natural logarithms 

Chapter 5. Factorials, Eulerian integrals, and related functions 

Chapter 6. Sine, cosine, exponential, and logarithmic integrals and related 
functions 

Chapter 7. Probability integrals and related functions 

Chapter 8. Elliptic integrals and elliptic functions 

Chapter 9. Legendre polynomials and Legendre functions 

Chapter 10. Cylindrical functions 

Chapter 11. Certain special functions and integrals 

Chapter 12. Solutions of certain equations 

Chapter 13. Sums and quantities connected with finite differences 

Chapter 14. Mathematical constants 

Chapter 15. Prime numbers, factors, products, quotients, and fractions. 


The user must be willing to thumb through several pages of the table of 
contents to find the table he seeks. This requirement might have been alleviated 
by an index of functions listed, such as the one in [1], and the reviewer would 
certainly be happier with such an index. However, this omission may amount to 
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a negligible time loss. Fresnel integrals, for example, are listed in Chapter 7. 
The user might seek them here in any one of several ways: adequate knowledge 
of special functions and the arrangement of the Guide might lead him to look here, 
he might leaf through the eighteen pages required before he finds them in the 
table of contents, he might look for Fresnel in the name index and find a reference, 
or he might seek in the name index the name of the author of some known table 
of Fresnel integrals. 

Bibliographic references are grouped together in the back of the book accord- 
ing to the chapter. Thus, under Fresnel in the index one finds a reference to the 
bibliography of Chapter 7. Under the chapter groupings, the bibliographic refer- 
ences are arranged alphabetically according to author or title (whichever seemed 
more appropriate to the authors) with Cyrillic characters preceding Latin. Page 
references are usually given to papers in journals, but not to collections of tables, 
such as Jahnke-Emde. Principal libraries in USSR owning the tables listed are 
noted in the bibliography. 

Typography in this book, at least so far as formulas and important listings 
are concerned, is excellent. A few letters are printed poorly in the verbal portions, 
but on the whole the book is easy to read, easy to use, and appealing. 

Utility of indexes of this kind must be obvious. 

Throughout the book notation is seemingly consistent with that of [1]. 
Perhaps one of the great services that [1 ] and this book may provide (along with 
the Bateman Project volumes [4 ]) will be a general consistency in notation among 
various users of special functions. It is unfortunate that this consistency has not 
developed spectacularly during the ten years that [1] has been available. The 
reviewer still hopes that the large efforts of the USSR table makers, the Royal 
Society table markers, the NBS table makers, and other large groups, will lead to 
some unification of terminology. At least the three gréups mentioned are largely 
consistent. However, we have recently seen Harvard University tables of the 
“Error Integrals’ and the Cambridge University Press edition of ‘‘Fresnel 
Integrals” listing functions which do not conform with the suggestions of these 
scholarly indexes. 

This Guide has an important place on the reviewer's bookshelf. 


i Ge 


1. A. FLetcHer, J. C. P. Miter, & L. RosENHEAD, An Index of Mathematical Tables, Scien- 

tific Computing Service Limited, London, 1946. 
- Kar Scuittte, Index Mathematischer Tafelwerke und Tabellen, R. Oldenbourg, Miinchen, 

1955. 

3. Harotp T. Davis & VERA FisHER, A Bibliography of Mathematical Tables, Northwestern 
University, Evanston, Illinois, 1949. 

4. ARTHUR ErpEéLYyI, WILHELM MaGnus, Fritz OBERHETTINGER, & FRANCESCO G. TRICOMI, 
Higher Transcendental Functions, McGraw-Hill Book Co., Inc., New York, 1953. 


50[F].—D. H. Lenmer, “On the diophantine equation x* + y? + 2 = 1,” 

London Math. Soc., Jn., v. 31, 1956, p. 275-280. 

The author makes a study of the diophantine equation (1), x* + y' + 2 = 1, 
and finds a sequence of explicit parametric solutions of which (2), x = 9%, 
y = — 94 + 3t, 2 = — 9 + 1, is the simplest. It is shown how some of these 
solutions give rise to an infinity of others. Several numerical instances are given, 
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of which the following are typical: 


x = 3753, y = —2676, = —3230, 
x = 3753, y = —5262, z = 4528, 
x = —98196 140, y = 78978 818, z = 76869 289. 


The study was prompted by a table of Miller and Woollett [1] in which of 21 
non-trivial solutions of (1) with |x|, |y|, |z| << 3164, only 8 are covered by (2). 
Unfortunately the solutions obtained by the author all fall outside the range of 
this table, and so there are still solutions of (1) unaccounted for by any parametric 
formula. 

Morris NEWMAN 
National Bureau of Standards 
Washington, D. C. 


J. C. P. Mmter & M. F. C. Woo.tertt, “Solutions of the Diophantine equation 
ot pt sek” London Math. Soc., Jn., v. 30, 1955, p. 101-110. 


51[F].—H. J. Gopwin, “Real quartic fields with small discriminant,’ London 

Math. Soc., Jn., v. 31, 1956, p. 478-485. 

The author lists all totally real quartic fields of discriminants less than 11,664 
giving full details (for simplicity) in the case of discriminants less than 2,000. 
There are 45 such fields with quadratic sub-field, starting with R(7 + 2-5#)! of 
discriminant 725, and 19 such fields without quadratic sub-field starting with 
one generated by x* — 8x* + 20x? — 17x + 3 of discriminant 1,957. The author 
conjectures the first two different fields of like discriminant are those generated 
by x* — 10x? + 30x? — 32x + 10 and R(17 + 4- 2+)! of discriminant 16,448. The 
tables corroborate B. N. Delone and D. K. Faddeev [1], who produced a similar 
table with discriminants up to 8,112. The computations were apparently done 
by hand. 

HarvEY COHN 
Washington University 
St. Louis, Missouri 


1. B. N. Detone & D. K. FappEev, “Teorifa irrafsional’nostei trel’ei stepeni” (Theory of 
cubic irrationalities), Akad. Nauk SSSR, Mat. Inst. Steklova, Trudy, v. 11, 1940. 


52[F, L].—D. H. Lenmer, “On the roots of the Riemann zeta-function,” Acta 

Mathematica, v. 95, 1956, p. 291-298. 

The author gives the results of numerical computations undertaken on SWAC 
relating to the behavior of the Riemann zeta-function {(s), s = ¢ + it, on the 
critical line ¢ = 1/2, ¢ > 0. The first 10,000 zeros of ¢(s) were calculated, and it 
was found that these are all, simple and have real part 1/2. Thus the Riemann 
hypothesis is true at least for ¢ < 9878.910. 

As an auxiliary tabulation, the number of failures of Gram’s law among the 
first 2 Gram intervals for » = 1000(1000)10,000 was recorded. Gram’s law fails 
for about 814 instances below 10,000. 

The author adds in proof that a few more hours of machine time were used 
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to examine the next 15,000 roots. The behavior of ¢(s) becomes steadily ‘‘worse’”’ 
but all roots have real part 1/2. 
The total computing time amounted only to several hours. 


Morris NEWMAN 
National Bureau of Standards 


Washington, D. C. 


53[H, X].—W. L. Witson, Jr., Tables of Inverses to Laplacian Operators over 

Triangular Grids, two typewritten leaves 83 X 11 inches and one transparency 

27 X 48 inches deposited in UMT FILE. 

Tables to 6D of the inverses of matrices L,,m = 2(1)7, where L, is a matrix 
of n(n + 1)/2 rows and columns which is an analogue of six times the Laplacian 
operator over a regular triangular grid of m(m + 1)/2 points. 

The Laplacian operator replaces a vector over the grid by itself minus one-sixth 
the sum of its components at the six or fewer adjacent points in the grid. 

These may be interesting in connection with [1 ]. 

W. L. WItson, Jr. 


University of California 
Los Angeles, California 


Calculation on SWAC was supported by the Office of Naval Research. 
1. I. J. Goon, “On the numerical solution of integral equations.’ [MTAC, this issue, p. 81-83.] 


54[I, X].—NBS Applied Mathematics Series, No. 35, Tables of Lagrangian 
Coefficients for Sexagesimal Interpolation, U. S. Govt. Printing Office, Wash- 
ington, D. C., 1954, ix + 157 p., 26 cm. Price $2.00. 

This work tabulates Lagrangian interpolation coefficients for use with tables 
with increments of one degree, one hour, or convenient fractions or multiples of 
these increments; entries are to 8D (with accuracy estimated at about 2 in the 
last decimal). In each case entries are for range 6 = 0(1’’)60’. 

These tables form a companion to the decimal tables [1 ], in which importance 
of such works is estimated. No comparable tables for sexagesimal arguments are 
known to the reviewer. 

Sixty randomly chosen arguments were tested in each table using SWAC 
(through the support of the Office of Naval Research), and the estimate of 
accuracy to within 2 in the last decimal is consistent with the SWAC results. 

Format is the standard style of the AMS, and the tables are, for the most part, 
easy to read and to use. Some faulty printing has been discovered, and a clarifica- 
tion sheet issued by the National Bureau of Standards for inclusion with the 
volumes lists the following values difficult to read in some copies: 


Page Minutes Seconds Coefficients 
23 43 8 A;: 0.6178 4506 
134 6 13 Ao: 0.9525 3516 
14 0.9523 7455 
15 0.9522 1377 
16 0.9520 5282. 


C.. B.:F. 


1. NYMTP, A. N. Lowan, technical director, Tables of Lagrangian Interpolation Coefficients, 
New York, Columbia University Press, 1944. [RMT 162, MTAC, v. 1, p. 314-315.] 
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55[1, X].—H. E. Sauzer, “Coefficients for complex osculatory interpolation over 
a Cartesian grid,” J. Math. Phys., v. 35, 1956, p. 152-163. 
The Hermite interpolation formula for a function f, given together with its 
derivative, at 21, ---, Zn, is 


f(z) = DL {Le (2) }2{1 — 2Le™’ (ee) (@ — 24) } fF (ee) + DO {Le (2) }? (2 — 24) f’ (ze) 
where the summation is for k = 1 to k = m and where 


Ly™ (2) = TT (2 — 2;)/ (ee — 25); 


j=1 


the error is a multiple of f°”. This is very handy when f and f’ are conveniently 
available, e.g., in tables where Jo and J; are given. 

In practice, we assume the z% are (some of) the vertices of the square in the grid 
containing z. In the case m = 2, we take 2 as the south-west vertex, 2; = 2 + h; 
if nm = 3 we add 22 = 20 + th and if nm = 4 we also take 2; = 29 + (1 + i)h. We 
have z = 29 + (p + ig)h where 0 < p; g < 1. The basic formula can be written 


f(%o + (b + ig)h) = LX {Ar (b + ig) f(z) + Bi (p + ig) f’ (ze) } 


where the A;™, B,™ are polynomials of degree 2m — 1 in p + ig. 

The tables give values of the polynomials A.™, B,™ for p = 0(.1)1, ¢ = 0(.1)1, 
for n = 2, 3, 4. For = 2, values are given to 3D; for m = 3, values are given 
to 5D; for m = 4, values are given to 8D; all values are exact. 

The construction and checking of the tables is described. 

5. a 


56[L ].—O. EMERSLEBEN, “‘Anwendungen zahlentheoretischer Abschatzungen bei 
numerischen Rechnungen,” Z. angew. Math. Mech., v. 33, 1953, p. 1-3. 
O. EMERSLEBEN, “Uber Summen Epsteinscher Zetafunktionen regel- 
massig Verteilter ‘unterer’ Parameter,’’ Math. Nachr., v. 13, 1955, p. 59-72. 
O. EMERSLEBEN, ‘Uber Zwei Epsteinsche Zetafunktionen 4. und 8. 
Ordnung,” Ber. Math. Tagung, 1953, p. 233-250. 
O. EMERSLEBEN, “Uber eine doppelperiodische Parallelstrémung zaher 
Fliissigkeiten,’’ Z. angew. Math. Mech., v. 35, 1955, p. 156-160. 
The author lists a large number of properties of the Epstein zeta-function in 
these papers, defined by 


(s) = 2 (—1)* 


n=1 


‘Tp (n) 
ni2 * 








0 
() 2|¢ 


where r,(m) is the number of representations of m as a sum of p squares, as well as 
some related functions. It is shown how a judicious use of the number-theoretic 
properties of the function can be used in numerical evaluation. Thus, the author 
gives to 6 decimals 


@2|3| (p), for p=1,2,3,4,6,8. 
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In addition, the author gives exact values or 6 or more decimal place values 
for the Epstein zeta-functions 

















0 0 0 0 lo 0 0 0| 0 0 
Zl $]@ 2]p 31 2[g g]@ 2] G|@ |p 4]@ 
0 0 0 0| 00 0 00 0 10 
Z\4 $| @. Zia 3] Zia a gf Zia a 3] 2/3), 

















as well as for some other functions which are combinations of these. 


Morris NEWMAN 
National Bureau of Standards 
Washington, D. C. 


57[L].—T. H. Soutuarp, “Approximation and table of the Weierstrass g func- 
tion in the equianharmonic case for real argument.’’ [MTAC, this issue, 

p. 99-100. ] 

A sparse table giving f(z; 0,1) = g(z;0,1) — 1/z2*, e = 0(.1).8(.05)1.55, 7D, 
with modified second differences. Interpolation accurate to 6D using Everett's 
formula is claimed possible. 

The paper also contains an approximation of H(y) = x* g(x;0,1) where 
y = x®§,0 < x < 1.53, 7S. This approximation is a cubic polynomial in y. 

Early tabulation to 5D is in [1] for the g function. 

[ pee 


1. EUGENE JAHNKE & Fritz Empe, Tables of Functions with Formulae and Curves, Dover 
Publications, New York, 1945, p. 102-104. 


58[L].—F. H. HoLLanper & C. B. Tompxins, “What you should know about 
digital computers,” Chem. Eng. Prog., v. 52, 1956, p. 451-454. 
T. H. Southard’s approximation (but not the table) to g(x; 0, 1). (See 
Review 57 above.) 
> a 


59(L, P, S].—G. Betrorp, L. Jackson Lastett, & J. N. SNyDER, Table Pertain- 
ing to Solutions of a Hill Equation, 1956, 420 p., two 8}” X 11” multilithed 
copies in notebook binders deposited in the UMT file. 

Solutions of a Hill equation of the form 


d*y/d? + (A + Bcos 2t + Ccos 4¢ + D cos 6t)Y = 0, 


have been tabulated for various values of the coefficients A, B, C, and D, as 
described by Belford, Laslett, and Snyder in this issue of MTAC [1]. 

The introduction to the table gives a derivation of the equations used to 
evaluate the tabulated quantities, and a detailed discussion of the method of 
estimating the truncation error. 


F. H. HOLLANDER 
University of California 
Los Angeles, California 


1. G. Betrorp, L. oe Lastett, & J. N. Snyper, “Table pertaining to solutions of a 
Hill equation.” [MTAC, this issue, p. 79-81. 
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60[M ].—D. Brerens DE Haan, Nouvelles Tables d’Intégrales Définies, Hafner 

Publishing Co., New York, 1957 (edition of 1867—corrected), xiv + 716 p., 

24 cm. Price $12.50. 

This is an unaltered reproduction of G. E. Stechert’s 1939 edition of the same 
tables, which in turn was a “corrected” facsimile edition of the original 1867 
edition of Bierens de Haan’s classical work, augmented by an English translation, 
due to J. F. Ritt, of the preface. No indication of the nature and extent of the 
corrections made can be found in either the 1939 or the present edition. 

Bierens de Haan’s gigantic work hardly needs recommendation. For the 
uninitiated reader we mention that it is probably the most extensive published 
collection of definite integrals involving elementary functions; there are no indefi- 
nite integrals, no definite integrals which can be derived from elementary in- 
definite integrals, and no integrals involving higher transcendental functions. 


PETER HENRICI 
University of California 


Los Angeles, California 


61[P].—Ernst GLowatTzkI, Sechsstellige Tafel der Cauer-Parameter, Abh. Bayer. 

Akad. Wiss. Math.-Nat. KI. (N.F.), v. 67, 1955, 37 p. 

These tables are of use in the design of electrical networks, see W. Cauer [1]. 
They give, for @ = 0(1°)90° and for m = 1(1)m, m = 1(1)12, the values to 
6D of a, = (sin 6)'sn (mK/n; sin 6), to 6S of A and to 3D of —InA. Here, 
if n is odd, A = a,"I1,a2,_1, where vy runs from 1 to $(m + 1), while if is even 
A = I,a*2,_1, where » runs from 1 to 4m. When n = 1, 2, 3, 5, 6, 9, 10, the compu- 
tations were carried out by hand to 12D, using the tables of G. W. and R. M. 
Spenceley [2]. When m = 4, 7, 8, 11, 12, the computations were carried out on 
the Géttingen automatic computer G1, using the series expansions, and working 
to 9D. No details of the checking are given. The tables are clearly printed. 

The parameters a,, are the zeros of the Zolotareff functions of degree , based 
on the interval J: (—&!, k#), where @ < k = sin @ < 1. The Zolotareff function is 
that rational function of x, with numerator and denominator of degree n, with 
value unity at x = 1, which gives best approximation (in the sense of minimal 
deviation) to zero in the interval J and to © in the complement of the interval 
(—k-+, k++). The maximal deviation from zero in J is A. For an account of these 
functions see H. Piloty [3]. 


p Soe 
This review was prepared by John Todd for Mathematical Reviews. 


a W. Caver, Theorie der linearen Wechselstromschaltungen, Akademie Verlagsgesellschaft, 
pig, ee 

2. G. W. & R. M. SpENCELEY, Smithsonian Elliptic Functions Tables, Smithsonian miscel- 
or collections, v. 109 (Publication 3863), The Smithsonian Institution, Washington, D. C., 


> H. Pmory, “‘Zolotareffsche rationale Funktionen,” Z. angew. Math. Mech., v. 34, 1954, 
p. 175-189. 


62[P].—T. Sasaki, The Effect of Earth and Sea on the Propagation of Radio Waves, 
Numerical Computation Bureau, Report No. 6, Tokyo, 1952. 14 mimeo- 
graphed pages of theory, 25 cm., and 6 tables of functions, 25 K 34 cm. 

The mathematical approach is similar to that of Sommerfeld [1]. The sea is 
treated as finitely conductive as well as the land. Diffraction of radio waves in 
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the air by the boundary of the land and the sea is calculated. The attenuation of 
radio waves in the air above the land and above the sea is estimated. The dipole 
is assumed to be vertical, and is placed at the surface of the land and of the sea. 
The results are plotted in terms of the electric field as a function of the distance 
from the land-sea boundary. Results are given for the case where the receiver is 
located on the surface of the earth and at higher altitudes. 


C. T. Tomizuka 
University of Chicago 
Chicago, Illinois 


1. A. SomMMERFELD, ‘“‘Uber die Ausbreitung der Wellen in der drahtlosen Telegraphie,”’ Ann. 
Physik, v. 28, 1909, p. 685. 


63(K, L, P].—J. Hatcoms LANInG, JR., & RicHarp H. Battin, Random Processes 
in Automatic Control, McGraw-Hill Series in Control Systems Engineering, 
McGraw-Hill Book Co., Inc., New York, 1956, ix + 434 p., 23 cm. Price 
$10.00. 

This book deals with the theory of random signals and noise as applied to 
automatic control system analysis and synthesis. 

Appendix B contains two tables of a functional relationship discussed earlier 
in the book. These two tables relate the ‘‘correlation’’ functions, respectively, of 
the input and output of two types of signal amplitude limiters. The input is 
supposed to be a stationary Gaussian stochastic process with zero mean and 
normalized correlation function, p(t). The limiters are supposed symmetric. The 
authors denote by x the ratio of the limiter threshold to the root mean square 
amplitude of the input. In the cases considered, for each value of ¢, the normalized 
output correlation functions depend only on the values of p and x. Table I gives 
the value of the normalized correlation function for the output of an “‘ideal” 
limiter for p = .05(.05)1.00 and x = .2(.2)2.0. Table II gives the corresponding 
relationship for the output of an ‘“approximate’”’ limiter for p = .1(.1)1.0, and 
x = .6(.2)2.0. 

Appendix D contains some properties of the polynomials in p~*t which are 
orthonormal functions with uniform weight on (0, ©). Tables are given of the 
unnormalized (i.e., multiplied by c~'/2) functional values of the first five such 
functions (i.e., those having degree 1, 2, 3, 4, and 5), for ct = 0(.01)1(.02)5(.1)9.9. 


HAROLD Davis 
University of California 
Los Angeles, California 


64[S].—J. L. Wotrson & H. S. GELLMAN, Five Figure Tables for the Conversion 
of Electron Momentum to Electron Kinetic Energy from 100 to 20000 Gauss-Cm, 
Atomic Energy of Canada, Ltd., Chalk River, Ontario, 1954 (reprinted July, 
1956), A.E.C.L. Report No. 327, i + 43 p., 27 cm. Price $1.00. 

The tables give the kinetic energy T in kev of an electron of momentum Bp 
in gauss-cm, according to the formula 


2 
T = 510.969 [V(=25) +1 - | 
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for Bp = 100(1)2000(10)10000(100)20000, 5S. Tables of proportional parts allow 
linear interpolation with an accuracy of 0.1 to 0.01%. 

Similar tables may also be found in Kai Sieghahn [1], and in NBS Applied 
Mathematics Series, No. 13 [2]. 


M. A. MELKANOFF 
University of California 


Los Angeles, California 


1. Kat SIEGHAHN, Beta- and Gamma-Ray Spectroscopy, North Holland Publishing Co., Amster- 
dam, 1955, p. 926-931. 

2. NBS Applied Mathematics Series, No. 13, Tables for the Analysis of Beta Spectra, U. S. 
Gov. + Office, Washington, D. C., 1952, p. 15-16. [RMT 1117, MTAC, v. 7, 1953, p. 
180-181. 


65(S, T ].—O. EMERSLEBEN, “‘Der Wert der Madelung Konstanten des Steinsalz- 
gitters,” Wissenschaftliche Zeitschrift der Universitat Greifswald, v. 3, 
1953-54, Math. naturw. Rethe, p. 607-617. 


This is a useful survey of the historical development of our knowledge of the 
accurate value of Madelung’s Constant (i.e., the constant A when writing the 
electrostatic energy per ion pair in an infinitely extended NaCl type lattice of 
alternating charges +e in the form Ae*/a, where a is the distance between nearest 
ions (=half the cubic lattice constant), and the ions are considered as point 
charges). The value first calculated by Madelung in 1918 is given in (1)—correct- 
ing a trivial numerical error—as 1.746(6). Ewald, employing a method of evaluat- 
ing the lattice sum }> + e*/r, by Fourier (or Theta-) transformation, obtained 


1.747 in 1921; and Emersleben, in considering the more general lattice sums 
x + e&/r,* as functions of the exponent s and by following in the complex s-plane 


the analytical continuation from the region of real s > 3, where convergence is 
absolute, to smaller real s values and eventually to s = 1, obtained in his disserta- 
tion with Born, 1922, a new approach to the mathematical convergence problem. 
The numerical value then given to 5 decimals was later (1950) extended by him 
to 14 decimal places, and independently by Y. Sakamoto (1953) in Japan to 16 
decimal places: A = 1.74756 45946 33182 163 + 9. The paper under review 
reveals very little about the details of the calculations, tables and machines used, 
or number of terms required, except indicating that the calculation followed 
essentially Ewald’s method and that special tables (with 15 and more decimal 
places) were calculated for the Error Integral with argument Vnz, n integer [1]. 

An interesting paragraph deals with the relation of the Madelung Constant 
and the Epstein Zeta-functions in p dimensions (see Review 56, p. 109), 


(p)Z|°| (s) = X ersery (nar, 


n=l 
where r,(m) is the number of lattice vectors of length Vm = (ny? + mn? +---+ ,?)! 
from the origin of a p-dimensional cubic lattice. It had been shown by Emersleben 


in 1950 that in one dimension — A = Z 1 (1) = 2 In 2 and in two dimensions 


-A=%Z : (1); furthermore that for p = 1, 2, 4, 8, it is possible to express 
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wz|3 (p) in a closed form (= — wIn2 for » = 2) and that there exist func- 


(s) for k = 4and k = 0 when p = 1, 2, 4, 8, 





tional equations between the (p)Z . 
such as 


(s) = (2!-#2 — 1). .@z\ 


7 | 0 
2), 9| (9): 














This last relation can conveniently be used for the calculation of the Madelung 
constant for a square net. These relations go back to relations between the num- 
bers r,(), i.e., the number of lattice points at distances ¥m from the origin. 
These relations are only known for p = 2, 4, 8, namely, if « denotes an odd value 
of n, and / any integer 


ro(2'u) = re(u); r4(2'U) = 3rg(u); r9(2'u) = 1/7(8.2%* — 15)rg(u). 


The search for further relations of this kind, which may lead to further relations 
between the Zetas and be used in calculating the Madelung constants is being 
continued by the author [2]. 


P. P. EwWaLp 
Brooklyn Polytechnic Institute 
wast New York 


. Orto EMERSLEBEN, ‘“‘Numerische Werte des Fehlerintegrals fiir -{nx,"" Z. angew. Math. 
J m7 v. 31, 1951, p. 393-394 
2. Orro 'EMERSLEBEN, “Ober Summen Epsteinscher Zetafunktionen regelmassig Verteilter 
‘unterer’ Parameter,” Math. Nachr., v. 13, 1955, p. 59-72. 


66[S, T].—O. EmersLEBEN, “Mit welcher Genauigkeit bendtigt man die Madel- 
ung-konstante eines Kristalls beispielweise vom Steinsalzgitter?’’ Z. fiir 
Physik. Chemie, v. 204, 1955, p. 43-55. 


In this paper the author stresses that the knowledge of an accurate value of 
Madelung’s constant A is of practical importance for a discussion of the Coulomb 
part of the surface energy of finite crystals, for the application of Born’s cycle 
(where all energy terms should add up to zero), and in determining which of 
several crystal structures is the more stable (e.g., NaCl or CsCl type). In these 
applications the uncertainty in the values of some of the other physical constants, 
such as electron charge, Avogadro number, or lattice constant, cancel away, so 
that it becomes important to know A correctly to about 10~ of its value. 

The paper contains a plot of Born’s ground potential which may be found 
useful for calculations of energies in lattices. 


‘ P. P. EwWALp 
Brooklyn Polytechnic Institute 


Brooklyn, New York 


67[L, X].—A. Erpéty1, W. Macnus, F. OBERHETTINGER, & F. G. TRICOMI, 
Higher Transcendental Functions. Based, in part, on notes left by H. BATEMAN. 
McGraw-Hill, New York, Toronto, London, v. 1, 1953, xxvi + 302 p., 23 cm. 
Price $6.50; v. 2, 1953, xvii + 396 p., 23 cm. Price $7.50; v. 3, 1955, xvii + 292 
p., 23 cm. Price $6.50. 

The late Harry Bateman left at his death in 1946 extremely voluminous 
materials for a planned Guide to the Functions of very wide scope. Consideration 
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was given to the possibility of using this material for the benefit of the mathe- 
matical world. The Bateman manuscript project has been supported by the 
Office of Naval Research and the California Institute of Technology, of which 
Bateman was such an eminent member. In the event, the four distinguished 
mathematicians principally engaged on the project, with Erdélyi as director, 
found it necessary, if a useful work was to be produced in a reasonable time, to 
narrow the scope of their efforts. They also, understandably, found it easier to 
compile their own account and to make only limited use of Bateman’s notes. 
What we have, therefore, is a joint work produced to the memory of Bateman, 
rather than an execution of Bateman’s design. The whole work consists of the 
three volumes now under review, together with two volumes of Tables of Integral 
Transforms reviewed separately [Rev. 107, MTAC, v. 10, 1956, p. 252-254]. 

It is clear that the volumes on Higher Transcendental Functions are of the 
first importance. Their general scope may be indicated by the chapter headings. 


Vol. Chap. 


I I The gamma function 
II The hypergeometric function 
III Legendre functions 
IV The generalized hypergeometric series 
V Further generalizations of the hypergeometric function 
VI Confluent hypergeometric functions 
VII Bessel functions 
VIII Functions of the parabolic cylinder and of the paraboloid of 
revolution 
IX The incomplete gamma functions and reiated functions 
X Orthogonal polynomials 
XI Spherical and hyperspherical harmonic polynomials 
XII Orthogonal polynomials in several variables 
XIII Elliptic functions and integrals 
XIV Automorphic functions 
XV_ Lamé functions 
XVI Mathieu functions, spherical and ellipsoidal wave functions 
XVII _ An introduction to the functions of number theory 
XVIII Miscellaneous functions 
XIX Generating functions 


A glance at this list is enough to show that the field covered, despite inten- 
tional limitation, is still unusually wide. As a handbook or compendium, in which 
proofs are normally omitted or indicated in outline, the work is distinguished not 
only by the breadth of its conception but also by the depth of its treatment. It 
discusses topics and contains results which are not touched upon at all in any 
other general compendium of comparable size. The references for each chapter 
are gathered together at the end of the chapter; completeness is rightly disclaimed- 
but the select references provided will be found highly useful. Any applied mathe, 
matician accustomed to work with Whittaker and Watson's Modern Analysis 
at hand may be encouraged to add Higher Transcendental Functions to his 
apparatus. 
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Prospective users may ask, Is this a completely reliable work of reference, in 
the sense that an applied mathematician may extract any formula from it and 
use it with confidence? There is no doubt that some members of the mathematical 
public expect an affirmative answer, but anyone with experience in such enter- 
prises is aware that the answer is likely to be negative, and so it is in this instance. 
It is no surprise to the reviewer to find a few errata given on loose sheets; he 
would have been astonished to find matters otherwise, and has no doubt that 
further errata will come to light. All past experience shows that the removal of 
all error from such copious and varied collections of results requires a prodigious 
effort, larger even than that which has already, to our great benefit, been lavished 
on the work. Trivial misprints are fairly easily found. Creditably reliable though 
the work is, the reviewer would nevertheless advise using it as a source of inspira- 
tion about the kind of results which exist, rather than as a collection of results 
guaranteed in every detail. 

It seems a pity that a work of such importance is not printed from type, but 
the varityping is good of its kind. 

The reviewer may perhaps claim licence to put on record that his first name 
is wrongly spelled (v. 2, p. 382, for Allan read Alan). 

A. F. 


68[X].—NBS Applied Mathematics Series, No. 15, Problems for the Numerical 

Analysis of the Future, U. S. Gov. Printing Office, Washington, D. C., 1951, 

iv + 21 p., 26.0 cm. Price $0.20. 

When asked to review a set of papers written over eight years ago on Problems 
for the Numerical Analysis of the Future, one is sorely tempted to seize the unfair 
advantage afforded by the passage of years and discyss just how far the writers’ 
prophecies have been fulfilled. Better feelings prevail, however, when the amount 
of work ia which one would be involved by such a display of ill-nature is realised; 
this reviewer is going to give a purely factual account of the contents of the 
pamphlet. 

The first paper, by D. R. Hartree, is entitled ‘Some unsolved problems in 
numerical analysis.”” Following a brief introduction designed to illustrate the 
difference between formal analytical and practical numerical methods, he gives 
examples of problems to which the answer was not known. These include the 
elimination of approximately known roots of polynomial equations, the solution 
of systems of simultaneous nonlinear algebraic equations, in which he considers 
the possibility of using a steepest descent method to minimise the sum of the 
squares of the residuals, and two problems concerning relaxation methods. Further 
sections concern characteristic value problems in ordinary differential equations, 
and the motion of a continuous distribution of charge under mutual forces be- 
tween its parts. Finally, he considers the psychological problems that may arise 
in numerical analysis and the use of computers. 

The contribution by S. Lefschetz is entitled ‘‘Numerical calculations in non- 
linear mechanics.’ He discusses what is known of the solutions of van der Pol’'s 
equation, and points out how valuable would be a numerical investigation for a 
whole range of values of u, the parameter occurring in the dissipative term. The 
necessity for numerical work is even more pronounced in the case of van der Pol’s 
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equation with a forcing term, in which the chief interest lies in the existence of 
subharmonic resonance oscillations. He ends with the hope that the advent of 
electronic machines will enable nonlinear problems to be tackled directly, without 
the customary process of linearization. 

Bernard Friedman, in ‘“‘Wave propagation in hydrodynamics and electro- 
dynamics,” considers examples of the propagation of waves in which the difficulty 
of obtaining a solution is due to (1) an inhomogeneous medium, (2) complicated 
boundary conditions, (3) nonlinear equations, and (4) free boundaries. The 
respective examples cited by him are (1) the propagation of electromagnetic waves 
through the atmosphere, (2) tides in the ocean, (3) flood waves in rivers, and (4) 
the impact of a sphere on a fluid. 

The fourth and final paper concerns the then relatively new subject of ‘‘Linear 
programing,” and is written by George B. Dantzig, a pioneer in this field. He 
discusses the problem in very general terms, illustrating it with references to the 
housewife buying food most economically, the Air Force requirements to use men 
and equipment in the most efficient manner, and the automobile industry's 
production problems. A brief description of the “‘simplex’’ technique is included, 
and he concludes with the observation that the computers then available were not 
fast enough to solve many of these problems with existing techniques. 


E. T. Goopwin 
National Physical Laboratory 


Teddington, Middlesex, England 


69[X].—ALBERT A. BENNETT, WitiiAM E. Mine, & Harry Bateman, Nu- 


merical Integration of Differential Equations, Dover Publications, Inc., New 

York, 1956, 108 p., 20 cm. Price $1.35. 

This is an unaltered reproduction of a monograph which appeared in 1931 as 
Bulletin 92 of the Committee of the Division of Physical Sciences of the National 
Research Council. Judging from the many references which have been made to 
the original edition, it has served well its purpose of easing the conscience of pure 
mathematicians who neglected the numerical aspects of differential equations. 
The book contains indeed much valuable information on topics which are more 
or less intimately connected with numerical computation; however, there is 
relatively little material on the subject of the title as it is understood today. 

Chapter I, ‘The interpolation polynomial” (40 p., including a bibliography 
of 13 p.), by A. Bennett, is an extremely scholarly account of the various classical 
representations of the Lagrangian interpolation polynomial, its integral and its 
first derivative. The subject of chapter II, “Successive approximations” (20 p.), 
by the same author, is given a much wider scope than would be necessary for the 
numerical solution of differential equations. There is a wealth of historical detail 
on such matters as the computation of ¥2 by Indian mathematicians, and the 
method of Taylor’s expansion is traced to its origins. From the practical point of 
of view the most useful part of the book is chapter III, “Step-by-step methods of 
integration” (17 p., including a bibliography of 3 p.), by W. E. Milne, a miniature 
version of the chapters 2, 4, 5, and 6 of Milne’s more recent textbook on the 
subject [1]. In chapter IV, ‘Methods for partial differential equations” (16 p.), 
by the late H. Bateman, difference methods are mentioned; naturally, nothing is 
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reported here on the problems which arise in their actual use (degree of con- 
vergence, stability, iterative methods for solving partial difference equations), 
because work in these lines was almost nonexistent in 1931. There are some re- 
marks, interspersed with many references, on variational methods, and more than 
a third of the chapter is taken up by a discussion of the extension of solutions of 
(ordinary) differential equations to non-integral values of a parameter, notably 
by E. T. Whittaker’s cardinal series. 
PETER HENRICI 


University of California 
Los Angeles, California 


1. W. E. MiLne, Numerical Solution of Differential Equations, John Wiley & Sons, Inc., 
New York, 1953. 


70[X ].—L. CourFiGnaL, Résolution Numérique des Systemes d’ Equations Linéaires, 

Gauthier-Villars, Paris, v. 2, 1956, iii + 180 p., 21 cm. Price $5.86. 

As the title implies, the author, in this volume, is concerned with the specific 
problem of finding the numerical solution to a system of linear algebraic equations. 

Chapter 1 contains the definitions and notation. The cracovians of Banachie- 
wicz [1] are used rather than the more familiar matrices. These tableaux use 
column-by-column multiplication instead of row-by-column multiplication and 
the author has reversed the usual meaning of the subscripts of an element a,j, i.e., 
the 7 denotes the column and the j denotes the row containing a,;. 

In the second chapter he introduces the concepts of the determinant of a 
tableau, the rank of a tableau, and the réduzt of a tableau. This latter term denotes 
a tableau obtained by a transformation which eliminates the mth row and mth 
column while replacing a;; by its transform, a’;;, where 


’ 
a’ ij = Aig — CinOmj/Omn- 


In the third chapter the author gives the formal solution of a system of linear 
equations. His algorithm involves the construction of réduits of a tableau, and the 
formulas are easily programmed for hand or machine computation. 

In chapter 4 numerical solutions are discussed, and a numerical analysis of the 
effects of errors in the coefficients and constants of the system is included. There 
is a discussion of round-off errors and methods to reduce these to a minimum. 
A Relaxation method of Southwell [2] is described briefly. The final section of 
the chapter is concerned with systems of equations which are not normal. 

The fifth and final chapter makes up approximately half of the book. It con- 
tains a detailed discussion of computational procedures which includes copies of 
computation forms useful for the hand computer and the computer using a desk 
calculator. Sample problems are worked out in detail, and the step-by-step pro- 
cedures used are presented in the form of lists of computation steps. These lists 
play the same role as the flow charts which some programmers construct prior to 
coding a problem for solution using an automatic computing machine. 

The author claims that the column-by-column multiplication employed using 
cracovians justifies his choice of notation. His viewpoint is undoubtedly influenced 
by the fact that the book is written as a manual for hand calculation. Automatic 
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calculation is mentioned only briefly (page 90), and the discussion is restricted to 
electronic analog machines. 

The book is well written and will serve as a good reference for those interested 
in manual methods for solving linear equations. The printing is easily legible and 
entirely satisfactory, although there are a few minor misprints, e.g., in the defini- 
tion of cracovian multiplication on page 12, the left member of equation (16) 
should be ||,¢|| rather than ||5,,||, and several page numbers listed in the table of 
contents are incorrect. 


RoBert T. GREGORY 
University of California 
Goleta, California 


1. T. BANACHIEWICz, “‘Résolution d’un systéme d’équations linéaires algébriques par division,”’ 
Enseignement Math., v. 39, 1942-1950, p. 34-45. 

2. R. V. SOUTHWELL, Relaxation Methods in Engineering Science, Clarendon Press, Oxford, 
1940. 


71{L, X].—J. Crank, The Mathematics of Diffusion, Oxford University Press, 

New York, 1956, vi + 347 p., 23.5 cm. Price $8.00. 

The author gives a very extensive and detailed discussion of solutions of the 
diffusion equation. The results are presented in closed form and graphically. The 
major portion of the book is concerned with analytical solutions. However one 
of the thirteen chapters is devoted to numerical methods for dealing with a 
parabolic partial differential equation. 

The first six chapters of the book are devoted to the derivation of the diffusion 
equation and its solution, in case the diffusion coefficient is constant. Separate 
chapters are devoted to the discussion of problems with plane, cylindrical and 
spherical symmetry, respectively. Both steady state and time dependent solutions 
satisfying various boundary and initial conditions are built up from the ele- 
mentary solutions of the diffusion equation. 

Chapter VII contains a discussion of a number of problems in which diffusion 
occurs in two distinct regions separated by a moving boundary or interface across 
which there is either a discontinuity in concentration or the gradient in con- 
centration, the dependent variable of the problem. The problems treated are 
variants of a single mathematical problem which was treated generally by P. V. 
Danckwerts [1 ]. 

Chapter VIII treats problems in which one substance may be absorbed by 
another through which it can diffuse and with which it can also react chemically. 
This chapter is concerned with two dependent functions of position and time, 
the concentration c, and the amount of immobilized solute s which are linked by 
a pair of partial differential equations one of which is a parabolic one. Methods 
for dealing with special equations of this sort are discussed. 

Chapter [X discusses solutions of nonlinear parabolic differential equations 
which may be reduced to ordinary differential equations by requiring the inde- 
pendent variables to occur only in particular combinations. Various methods of 
successive approximation to the resulting nonlinear ordinary differential equation 
are given. 

Chapter X, entitled Finite-Difference Methods, is concerned with methods 
for obtaining numerical solutions to parabolic partial differential equations. The 
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author discusses the equation 
ac/aT = dc/ax? 


as an illustrative example and in terms of it gives various “recipes” for obtaining 
numerical solutions. Although the critical ratio 5T/ (6X)? where 6T and 6X repre- 
sent the mesh sizes used in T and X respectively is restricted in each of the recipes 
to be less than one the author never describes adequately the role that this quan- 
tity plays in the convergence of the numerical solution to the analytic one and 
the stability of the numerical algorithm. Moreover the possibility of obtaining 
stable and convergent algorithms with 57/(6X)? > 1 in cases where suitable 
boundary conditions are given is not even mentioned in spite of the fact that 
implicit algorithms are given. 

This inadequacy in the discussion of numerical methods may be related to the 
fact that the author’s discussion of the direct use of high-speed digital machines 
in the solution of parabolic differential equations is restricted to two sentences 
although a page is devoted to King’s work on the application of the Monte Carlo 
method to diffusion problems and five pages are devoted to a discussion of the 
use of differential analyzers and other analog devices for the obtaining of 
numerical solutions to diffusion problems. 

Thus while this chapter contains some useful information concerning the 
obtaining of numerical solutions to diffusion problems, it cannot be said to be a 
modern or definitive discussion of the numerical analysis involved or even of the 
computing instruments available. 

The coatents of Chapters XI and XII are described by their titles which are 
The Definition and Measurement of Diffusion Coefficients and Some Calculated 
Results for Variable Diffusion Coefficients respectively. 

Chapter XIII is concerned with the problem of the diffusion of one substance 
through the pores of a solid body which may absorb and immobilize some of the 
diffusing substance with the evolution or absorption of heat. The mathematical 
problem involved is one in which there are two dependent variables C, the con- 
centration, and 7, the temperature, which satisfy a pair of coupled parabolic 
partial differential equations. 

The book contains twenty short numerical tables of various mathematical 
functions which are needed in the evaluation of expressions given in the text. 
Included in these is a table of the error function and associated functions and 
tables giving the roots of various transcendental equations. The methods of 
computing these tables are not discussed. 

A. H. T. 


1. P. V. DANcKWweERrts, “Unsteady-state diffusion or heat-conduction with moving boundary,” 
Faraday Soc., Trans., v. 46, 1950, p. 701-712. 


72[M, S, X].—Marcet vAN LarETHEM, Une Méthode Nouvelle et Générale de 
Calcul des Intégrales Généralisiées, Publications Universitaires de Louvain, 
Louvain (Belgique), 1956, viii + 180 p., 23 cm. Price 250 Fr.b. ($5.00). 
This book is concerned with defining, evaluating the integral, 
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where y(x) has a singularity, a pole, at x = b and applying the results of the 
evaluation to various problems. The author considers the quantity 


x) = [yas 


and evaluates x as a power series in x. This series is then inverted to give x as a 
function of x. The coefficients of the various power series are given in detail as are 
estimates for the remainders. 

The main application is to the problem of determining the trajectories of ions 
in a cylindrical condensor. The results of the method used by the author for the 
evaluation of integrals of the form of J, are compared to the results obtained by 
a method used by O. Godart which involves a change of variable and numerical 
evaluation of the resulting integral. 

A. H. T. 


73[Z |.—The Proceedings of the Institution of Electrical Engineers, v. 103, Part B, 
Supplement No. 1-3, 1956 (Convention on Digital Computer Techniques), 
British Institution of Electrical Engineers, London. Price one pound. 


A well organized picture of the state of digital computer techniques in Great 
Britain was presented in London during April 1956. The full proceedings of the 
“Convention of Digital Computer Techniques’”’ has now been published by the 
British Institution of Electrical Engineers, Savoy Place, London W C 2, as Part 
B, Supplement Numbers 1, 2, 3 of its Proceedings. 

The portion of the conference concerned with applications of digital computers 
was particularly strong in its examples of engineering applications. Eight papers 
described the methods and results of machine solutions of problems in proton 
synchrotron design, load flow in power distribution systems, power systems 
engineering switching circuits, transformer design, synchronons motor calcula- 
tions, electric traction, the design of linear and nonlinear control systems, and 
the solution of a nonlinear heat conduction problem concerned with the freezing 
of fish. In most cases the authors were cognizant of the tutorial value of these 
papers beyond the time of presentation and went to greater than average pains 
to describe the preparation of the problems for the digital computer. 

A relative scarcity of papers on business and industrial control applications 
was compensated for by extensive discussion, particularly concerning the indus- 
trial control area. 

Formal problems in numerical analysis were dealt with by a few papers con- 
cerned with the solution of linear elliptic partial differential equations, calculation 
of characteristic roots and vectors of a real symmetric matrix, predictor-corrector 
methods for numerical integration, and the general programming strategy de- 
veloped at the University of Manchester Mark I computing facility. 

The description of modern digital computer systems was separated into cate- 
gories of ‘Commercially Available Computers” and ‘Experimental Computers.” 

The commercially available machines presented were the ““DEucr”’ (English 
Electric Co.); ‘‘MERcuRY” and “PrGasus” (Ferranti); “‘400’’ Series (Elliott 
Bros.) ; ‘‘HEC” and a punched card calculator series (British Tabulating Machine 
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Co., Ltd.); PCC (Powers-Samas); and the ‘‘650,”’ ‘‘704,” and ‘‘705’’ (IBM). Of 
particular interest are the characteristics of the “‘DEuUCE’”’ and the ‘“‘MERcurRy.” 

The ‘‘DEUCE”’ is based on the development of the ‘‘PrLot AcE”’ at the Nationa! 
Physical Laboratory. Confidence derived from the high reliability and produc- 
tivity achieved with the prototype machine led to an increased level of speed and 
complexity of organization using the same acoustic delay line and circuit com- 
ponents. A graded sequence of acoustic delay lines containing 1, 2, 4, or 32 words 
are incorporated with a highway busbar system to provide very flexible organiza- 
tion of information transfers between the delay lines. These fast access memories 
are backed up by 8192 words of magnetic drum storage. Optimal programming 
techniques are facilitated by the flexible organization and instruction code 
structure. 

The “MeErcurRY”’ is based on the development of the University of Manchester 
Mark II. ‘““MERcuRY’”’ is a serial machine using a digit repetition rate of 1 Mc/s. 
A random access magnetic core memory is provided compatibly by parallel 
transfers of 10 bits with a 10 microsecond cycle. The random access memory is, 
in turn, backed up by four magnetic drums each having a capacity of 163,840 
bits. Seven B-Registers and efficient floating point and fast multiplication or- 
ganization lead to high operating speeds and a powerful instruction code. The 
description of ‘‘MERCURY’”’ restricts itself primarily to logical differences between 
it and the Mark II. Full design details are covered in other papers describing the 
experimental Manchester machine. 

The experimental machines discussed were MArK II (University of Man- 
chester); an acoustic delay line accounting machine (British Tabulating Machine 
Co.); ‘‘NicHoLas” (Elliott Bros.); ‘‘Epsac II’’ (Cambridge University) ; ‘‘Imp” 
(Imperial College, University of London); ‘‘Ace”’ (National Physical Labora- 
tory); ““-BESM” (USSR Academy of Sciences). As indicated earlier the Man- 
chester group has done a substantial amount of new work and provides an ex- 
cellent 21-page exposition of the Mark II organization. A concise description of 
the Epsac II organization gives a little more insight into the microprogramming 
approach of the Cambridge group. 

In the same category of experimental machines, there were two papers pre- 
sented in a session on “‘logical design.’’ The first of these provides an extensive 
discussion of the thinking leading to the organizational form of the projected 
“ICCE II” at Imperial College. This is a very valuable paper insofar as it de- 
scribes approaches often left unmentioned in expositions occurring after the fact 
of machine construction. The large number of variables and inconclusiveness of 
decisions reached however serve to emphasize the need for a studied characteristic 
set of problems which might be used as a coarse yardstick for relative evaluation 
of machine organizations. 

The second paper discussed an experimental machine called ‘‘TACT”’ being 
developed at the National Physical Laboratory. The main feature of this machine 
is a step in the direction of automatic optimal programming for a machine with a 
large main cyclic memory. A small random access memory is provided and num- 
bers are stored with a “‘tag,”’ i.e., an additional set of bits carried with the number. 
When a three-address instruction is to be performed, it first inquires whether any 
of the required numbers are in the immediate access memory cells. If so, the 
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operation may be performed immediately; if not, it transfers one of the immediate 
access quantities to the lower speed memory and brings the necessary quantity 
to the immediate access cell. In computations with high repetition rate of opera- 
tion on a few quantities this type of organization would lead to increased efficiency. 
Most of the latter half of the convention concerned itself with component 
researches. 
The session on rapid access storage described the magnetic core memory de- 


veloped for Ferranti’s ‘‘MERCURY,” magnetic core logical circuitry, multistable 
vacuum tube circuits, and the design of the CRT memory for The Manchester 
Mark II. 

In addition to two papers describing magnetic tape storage devices at Cam- 
bridge and Ferranti, there was a presentation of work on variable reluctance 
reading heads being done at Manchester University. The latter is an important 
area of investigation since it provides one avenue towards the alleviation of the 
bottleneck between the high-speed computer and terminal data presentation. 

A number of papers on the transistor as a computing element describe the 
design of and experience with two small computers at the Atomic Energy Re- 
search Establishment and at Manchester University. In addition to these there 
are two good papers presenting studies of circuit techniques using transistors and 
magnetic cores at Manchester. 

Two digital analog conversion schemes were discussed. The first utilized a 
binary number to time interval to voltage output conversion achieving an accu- 
racy of one part in two inches. The second used a binary coded disc to digitalize 
a rotating shaft position to an accuracy of the order of one minute of arc. 

Several points of special interest emerged in the session on the computer in a 
nonarithmetic role including game playing, pattern recognition, and mechanical 
translation. A. D. Booth mentioned a promising avenue towards mechanical 
recognition evaluated at Birkbeck College. A voltage time picture of sound is 
sampled and the number representing the waveform, its time derivative, and its 
time integral are established as the criteria for recognition. Booth reports an 
accuracy of better than 0.1% in recognizing the spoken digits 0 to 9. For character 
recognition he generates numbers representing the first intersection of the 
character with each raster line. A predetermined set of numbers for the characters 
to be recognized is then compared with the output of the circuit. A minimum of 
the sum of the absolute values of the differences is sought. If there is ambiguity 
in the result for a chosen discrimination interval, the system goes back and 
generates a second set of numbers based on another criterion and compares with a 
second precalculated set to remove the ambiguity. No quantitative results were 
reported. Another paper describes work at IBM in the U. S. on character recogni- 
tion, where the problem is attacked by providing a digital computer with the 
quantized representation of the entire character field. With this set of information 
the analysts are then in a position to test various scanning and testing procedures 
using the logic of the machine. Such studies led to the development of the par- 
ticular system described in which the character is scanned serially with a fine 
spot; the light signals are quantized into black and white cells; the binary data 
from each vertical scan line are coded in three decimal digits representing the 
number, position, and size of the black areas. A large distribution of character 
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samples was taken from inked ribbons used for 150% of nominal life with identical 
samples taken at equal intervals during that period. The criteria for recognition 
were then established by a “‘learning’”’ process of testing and modifying the criteria. 

Concerning mechanical translation, A. F. Parker-Rhodes of the Cambridge 
Language Research Unit presented the type of computer operations needed for 
further advances and I. S. Muchin of the Academy of Sciences of the USSR 
described in detail an experiment in machine translation carried out on BESM. 

The last session of the convention covered miscellaneous circuit techniques 
and components of computers. These included several papers on the use of electro- 
magnetic and acoustic delay lines. One interesting paper from the Manchester 
University group described a multi-input analog adder for use in a fast binary 
multiplier. Another paper provided a good exposition of the design of a parallel 
arithmetic unit at Imperial College. 

As a whole, the publication of the convention proceedings represents a contri- 
bution to the digital computer field not only by providing the results of British 
ingenuity, but by the effort made in many papers to describe their work without 
assuming an intimate knowledge of it by the reader. This effort strongly affects 
progress. 


GERALD ESTRIN 
University of California 


Los Angeles, California 


74[Z].—S. A. LeBepev, Elektronnye Vychislitel’nye Mashiny i Obrabotka In- 
formaisit (Electronic Computing Machines and Information Processes), Akad. 


Nauk SSSR, Moscow, 1956, 47 p., 20 cm. Price 65 kopeks. 


This pamphlet is one of the series of educational booklets issued by the 
Academy for popular consumption. Its author is the Academician in charge of 
BESM (this name is made up of the initials of the four Russian words standing 
for High-speed Electronic Computing Machine) and principal compiler of the 
Academy’s Index of Mathematical Tables, Moscow, 1956 (Rev. 49, p. 104-106, 
this issue). The booklet is couched in delightfully lucid language and is profusely 
illustrated with diagrams, photographs, and sample programs. The modest price, 
equivalent to about 15 cents, makes it available to anyone having an interest 
in the subject, and its high quality guarantees complete satisfaction to each 
purchaser. 

The Introduction points out the benefits of high-speed computation and 
emphasizes, in particular, the achievements made possible through the use of 
BESM not only in numerical computation but also in the initial—and quite 
spectacular—attempts to translate scientific English into Russian. 

The first chapter is devoted to a careful explanation of how electronic machines 
function. In the second, we find an excellent discussion of the binary versus the 
decimal systems with a complete explanation of the engineering principles which 
are used to represent the former system. This is followed by some very timely 
observations on the pros and cons of floating radix point representation of 
numbers. 

The last Chapter deals with a detailed description of BESM. Since the features 
of this machine are already familiar to the readers of this magazine, there is no 
need to repeat the facts here [1 ]. 
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The pamphlet is enriched by an excellent appendix, which can easily serve 
as a textbook for tyros in coding. It outlines the basic principles of program- 
preparation, and illustrates them by actual examples. It touches upon the use 
of subroutines and the choice of suitable checks to insure the accuracy of the 
computations. 

An altogether charming booklet for 65 kopecks! 


IpaA RHODES 
National Bureau of Standards 
Washington, D. C. 


1. S. A. LEBEDEV, “The high-speed electronic calculating machine of the Academy of Sciences 
of the U.S.S.R.,” J. of the Assn. for Computing Machinery, v. 3, 1956, p. 129-133. 


TABLE ERRATA 


252. A. A. Abramov, Tabliisy In [2] » Kompleknot Oblasti (Tables of in T[z] in 

a Complex Region), Akad. Nauk SSSR, Moscow, 1953. 

Published in the back of A. D. Smirnov’s book, Tabliisy funkisy Etri i 
Speisial'nykh Vyrothdennykh Gipergeometricheskikh Funkisy (Tables of Airy Func- 
tions and Special Degeneration of Hypergeometric Functions) (Akad. Nauk SSSR, 
Moscow, 1955) are lists of errors discovered in the tables of the Akad. Nauk series 
issued previously. Included in these lists of errors are the following in Abramov’s 
tables: 


As printed 
6 from right — 060 — 160 
3 from left 302 102 
6 from right 89 79 
3 from left 76 70 


Page Line Column 
14 1 from top 
26 18 from bottom 
54 19 from top 
61 17 from bottom 
82 24 from top 


Should be 


2 from left 397633 297633 


82 25 from top 


104 
131 
141 
143 
190 


11 from top 
8 from top 
10 from top 
1 from top 
4 from top 


2 from left 
4 from right 
5 from left 
4 from left 
1 from right 
3 from left 


202955 
49 

22 
673267 
80,12817 
29 


302955 
46 

32 
073267 
0,128178 
19 


217 5 from bottom 
218 2 from top 
273 15 from bottom 
286 8 from bottom 
319 17 from top 
321 8 from top 


2 from right 22 32 
2 from right 0 10 
2 from right 37 27 
5 from left 0 9 
4 from right 45 25 
6 from right 35 25 


253. AKADEMIIA Nauk SSSR. Institut Tochnol mekhaniki i vychislitel’noi 
tekhniki. Matematicheskie Tablitsy. Desiatiéna Tabliisy logarifmov kompleks- 
nykh chisel i perekhoda ot dekartovykh koordinat k poliarnym. Tabliisy funkisit 
[Ten place tables of logarithms of complex numbers and of the transformation from 
cartesian to polar coordinates. Tables of functions ] In x, arctg x, $ In (1 + x*), 
(1 + x*)#. Moscow, 1952. [See RMT 1206, MTAC, v. 8, 1954, p. 149.] 
Published in the back of A. D. Smirnov’s book, Tabliisy funkisy Eiri i 

Speisial’nykh Vyrozhdennykh Gipergeometricheskikh Funkisy (Tables of Airy Func- 
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tions and Special Degeneration of Hypergeometric Functions) (Akad. Nauk SSSR, 
Moscow, 1955), are lists of errors discovered in the tables of the Akad. Nauk 


series issued previously. Included in these lists of errors are the following for 
these tables: 


Page 


13 
16 
22 
23 
35 
35 
38 
56 
82 
98 
100 
101 
102 
106 
107 


Line 
10 from top 
6 from bottom 
25 from bottom 
8 from bottom 
8 from top 
6 from bottom 
23 from top 
22 from bottom 
8 from top 
8 from top 
24 from bottom 
15 from top 
13 from bottom 
9 from top 
18 from top 


Column 


2 from left 
2 from left 
2 from right 
2 from right 
4 from right 
2 from right 
1 from leit 
1 from left 
2 from right 
4 from left 
5 from right 
4 from left 
4 from left 
5 from right 
4 from left 


As printed 


6,4756128 680 
9,6647477 061 
2882 741 
8714 020 
4,857 

3567 723 
5,122 

9,928 

1118 566 
9890 

3330 

8729 

7822 

6267 

4595 


Should be 


0,4756128 680 
0,6647477 060 
3882 741 
3714 020 
3,857 

2567 723 
4,122 


Insert 6 from top 11 from left 0276 


254. AKADEMIIA Nauk SSSR. [Institut Tochnoit mekhaniki i vychislitel’nol 
tekhniki. Matematicheskie Tablitsy. Tablitsy Integralov Freneha (Tables of 
Fresnel Integrals), Moscow, 1953. [See Rev. 40, MTAC, v. 9, 1955, p. 75-76.] 
Published in the back of A. D. Smirnov’s book, Tabliisy funkisy Eiri i 

Speisial'nykh Vyrozhdennykh Gipergeometricheskikh Funkisy (Tables of Airy Func- 

tions and Special Degeneration of Hypergeometric Funetions) (Akad. Nauk SSSR, 

Moscow, 1955) are lists of errors discovered in the tables of the Akad. Nauk 

series issued previously. Included in these lists of errors are the following for these 

tables of Fresnel Integrals: 


Should be 


4857421 
4857681 
5010046 
5081338 


Column 


2 from right 
2 from left 
2 from left 
2 from left 


Page Line 

144 8 from top 
155 4 from top 
169 8 from bottom 
232 23 from top 


As printed 
4857424 
4857684 
5010045 
5081336 


255. ISABELLE ARSHAM, Chebyshev Coefficients for Chebyshev Polynomials of 
Orders 12 and 24 under the General Linear Transformation (U), Report No. 
TR-326, Diamond Ordnance Fuze Laboratories, Ordfiiance Corps, Department 
of the Army, Washington, D. C., 1956 [Review 72, MTAC, v. 10, 1956, p. 231]. 
The following errata have been discovered : 


Page 3, I. Introduction: formula (2) 


aa Ta(x) = "4 (in arc cos x), 


should read 


T(x) = cos (m arc cos x), 
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Page 4, III. Method of Computation: formula (3) 


[(n—m)/2] [(n—m)/2]—r 
reads bn = > > Cm--2j, n—m—2r—2j0"* 71h i 
r=0 i—0 
[(n—m)/2] [(n—m)/2]—r 
should read a= > DCm 24, 2 —2r—2 Ih 22, 


T=0 j=o0 


Page 4, line —6, 


$ 
reads where ¢ is the general - - - 


should read where (:) is the general - - -. 


Page 15, column 3, sixth entry, 


reads 28555887360a""b* 
should read 2855588736005‘. 


Note from DIAMOND ORDNANCE FuZzE LABORATORIES 


NOTES 
John von Neumann 


1903-1957 


John von Neumann died in Washington, D. C. on 8 February. It is probable 
that no other person had done more to advance the development and application 
of modern electronic computers to their present state. 

His contribution to this advancement was on at least three fronts; he realized 
the computational potential of components build of modern electronic com- 
ponents; he formulated and helped solve the real logical difficulties which stood 
in the way of exploiting this potential; he recognized and fearlessly attacked prob- 
lems of importance whose solution involved millions and tens of millions of 
multiplications. 

In another way his geniality, his reputation, and his most remarkable talent 
for clear exposition initiated, more than anything else, the acceptance, first by 
the Government and then by industry, of our present large program of construc- 
tion of machines. 

His interest in electronic calculators continued to the end of his scientific 
career. He expounded doctrines of increasing the size and complexity of these 
machines when componentry was available; also he studied the analogy between 
the operations they carried out and the operations of human thought. He con- 
sidered deeply the question of reliable computation with unreliable components. 
He contributed in countless other ways to the development of the computing 
devices presently available. 
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His interest was not only in the device, but in the applications, and many 
examples attest to his use of the new machines in deep and scholarly researches. 
In addition, he pioneered in the researches in numerical weather prediction, in 
Monte Carlo techniques, and in many other numerical procedures which are 
currently subjects of active research. So penetrating was his thought on these 
matters and so lucid his explanations that the present support of many of these 
research projects is almost wholly attributable to him. Much of the inspiration of 
the people continuing to work on these projects must also have been induced by 
von Neumann. 

Von Neumann’s remarkable expository talent seemed to lie in his ability to 
grasp and to state the most elementary principles upon which a problem depends 
and upon his ability to take an elegant and direct path from such statement of 
principles to the solution of the problem. I have never known of any problem 
attacked by him without an outcome of direct constructive suggestions; the field 
of the problem and the competence of the investigator approaching him seemed 
immaterial, for von Neumann still made valuable and profound contributions to 
all problems he agreed to examine. His speed in this was astounding. 

He was a man devoted to his ideals, and his last years were spent almost 
wholly in government service. At the time of his death he was a member of the 
Atomic Energy Commission, and he performed services for this Commission after 
his illness had advanced to a degree which made such service seem impossible. 
At the same time, so long as he was able, he continued his genial and illuminating 
advice on a profusion of topics covering almost the whole of mathematics and 
mathematical physics. 

Few men have lived who could apply themselves to as wide a range of prob- 
lems as von Neumann could; he seemed to bring immense power to bear instantly 
on any problem which interested him, and there was never any inertial delay in 
adjusting himself to the analysis of the new problem. 

His charm and graciousness were such that almost every person who knew 
him hoped to be his friend and was made to feel that he was a friend. This was 
partly due to the alert mind which always found something to analyze in a most 
fascinating way—he was interested in everything from abstruse mathematics and 
physics to our illogical social mores and our even less logical spelling. I remember 
him explaining to my very young son who had a large bump on his head that 
“Everybody falls on his head sooner or later, some just seem to have fallen a little 
too hard.” During the same visit he told my wife his theory of entertaining babies: 
“Just say ‘goo, goo’ and poke them gently here.’’ When the son was five he came 
on a conversation on our porch in which von Neumann was analyzing words 
which had some quaint touch when spelled backwards. Von Neumann was per- 
plexed when the son asked what OPPO spelled and he was fascinated with the 
answer—OPPO spells POOP inside out. My son, now seven, feels a real sense of 
loss in the death of von Neumann. 

It must go without saying that those of us who realize that we owe so much 
of our knowledge, our advancement, and possibly even our civilization to him 
feel an unspeakable loss. 

C. Bi. 





CORRIGENDA 


NBS-NSF Training Program 


The National Bureau of Standards has announced that a Training Program 
in Numerical Analysis for Senior University Staff Members is now under way in 
Washington, D. C., with support from the National Science Foundation. The 
dates are February 11 to June 9, 1957. The program is under the general direction 
of Mr. John Todd, Chief of the Numerical Analysis Section, Applied Mathematics 
Division, National Bureau of Standards, with the assistance of Dr. F. L. Alt, 
Dr. H. A. Antosiewicz, Dr. P. Davis, Dr. Ky Fan, Dr. M. Marcus, Dr. M. 
Newman, Dr. F. Oberhettinger, Dr. O. Taussky, and others. 

The participants, limited to twelve in number, are selected regular university 
staff members who will receive training to enable them to direct the operation 
of university computing centers, and to organize training and research in nu- 
merical analysis in their own institutions upon completion of this course. 

The training program has been divided into three parts: 

a) Two weeks—Introduction to numerical analysis and to programming for 
automatic computers, 

b) Thirteen weeks—Survey of various chapters in numerical analysis and 
observation of operation of a general purpose computation laboratory, 

c) One week—Administrative problems. 

The program consists of formal instruction by mathematicians for mathema- 
ticians. Special topics will be covered by guest lectures. 


From NBS Announcement 


Wayne State University 
Conference on Matrix Computations 


The Department of Mathematics of Wayne State University announces that 
a Conference on Matrix Computations will be held on campus on September 3-6, 
1957. The Conference will be concerned with the basic problems of matrix algebra, 
the mathematical methods used to solve systems of linear equations, to compute 
the inverse of a matrix, and to find characteristic values and characteristic vectors. 
Papers suggesting new methods for the solution of standard problems will be 
presented and new problems demanding the efforts of mathematicians will be 
stressed. Smaller groups with well defined common interests will form discussion 
panels in the afternoons. 

Further announcements may be procured from Professor Wallace Givens, 
Department of Mathematics, Wayne University, Detroit 2, Michigan. 


CORRIGENDA 


HARVARD UNIVERSITY, Computation Laboratory, Annals, Tables of the Func- 
tion arc sin z, MTAC, v. 10, 1956, Review 66, p. 229, line 3 


for Marcell read Marcel 
for d'Etudes read d'Etude; 
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page 229, line 6, 
for 1949 read 1948. 


Mark Lorktn, “A set of test matrices,” MTAC, v. 9, 1955, p. 161, lis 
from bottom, reference 8, 


for 1940 read 1950. 


M. ScHuLeR & H. GEBELEIN, (a) Acht- und neunstellige Tabellen zu 
elliptischen Funktionen (Eight and Nine Place Tables of Elliptical Functie 
(b) Féinfstellige Tabellen zu den elliptischen Funktionen (Five Place Tables 
Elliptical Functions), MTAC, v. 10, 1956, Review 102, p. 247, line 19 from botta 


for @ and 2 again read @ and g again. 








